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1.  Research  Project  Overview 

This  report  describes  the  progress  and  results  of  the  University 
of  Southern  California  image  processing  research  study  for  the  period  of 
1  March  1972  to  31  August  1972. 

The  image  processing  research  study  has  been  subdivided  into 
five  projects: 

Image  analysis  projects 
Image  coding  projects 

Image  restoration  and  enhancement  projects 
Image  detection  and  measurement  projects 
Image  processing  support  projects 

The  image  analysis  project  comprises  the  background  research  effort 
into  the  basic  structure  of  images  in  order  to  develop  meaningful  quanti¬ 
tative  characterizations  of  an  image.  In  image  coding  the  orientation  of 
the  research  is  toward  the  development  of  digital  image  coding  systems 
that  represent  monochrome  and  color  images  with  a  minimal  number  of 
code  bits.  Image  restoration  is  the  task  of  improving  the  fidelity  of  an 
image  in  the  sense  of  compensating  for  image  degradations.  In  image 
enhancement,  picture  manipulation  processes  are  performed  to  provide 
a  more  subjectively  pleasing  image  or  to  convert  the  image  to  a  form 
more  amenable  to  human  or  machine  analysis.  The  objectives  of  the 
image  detection  and  measurement  projects  are  the  registration  of  images, 
detection  of  objects  within  pictures  and  measurements  of  image  features. 
Tinally,  the  image  support  projects  include  research  on  image  processing 
computer  languages  and  the  development  of  experimental  equipment  for 
the  sensing,  processing,  and  display  of  images. 

The  next  section  of  this  report  summarizes  some  of  the  research 
project  activities  during  the  past  six  months.  Sections  3  to  6  describe 
the  research  effort  on  the  projects  listed  above  during  the  reporting  period. 
Section  7  contains  a  short  description  of  new  projects  that  are  being  ini¬ 
tiated,  and  are  not  yet  to  the  reporting  stage.  Section  8  is  a  list  of  publi¬ 
cations  by  project  members. 


2.  Research  Project  Activities 

The  following  sections  describe  some  of  the  significant  pro¬ 
ject  activities  of  the  past  six  months. 

Publications .  July,  1972  has  been  a  banner  month  for  publications 
in  the  field  of  image  processing.  The  Institute  of  Electrical  and  Elec¬ 
tronic  Engineers  published  a  special  issue  of  the  Proceedings  of  the 
IEEE  on  digital  image  processing.  Professor  Harry  C.  Andrews  of 
use  acted  as  co-editor  of  the  special  issue.  The  journal  contained 
five  papers  written  by  the  USC  staff.  The  IEEE  also  devoted  the 
July  issue  of  the  Transactions  on  Computers  to  the  subject  of  two- 
dimensional  digital  signal  processing.  USC  staff  members  contri¬ 
buted  two  papers  to  the  issue.  Finally,  the  IEEE  carried  a  cover 
article  on  digital  image  processing  featuring  the  USC  research  on 
pseudocolor  image  enhancement.  A  reprint  of  this  paper  is  in¬ 
cluded  in  Section  4.  7. 

Image  Coding  Conference.  The  University  of  Southern  California 
has  been  chosen  as  the  host  for  the  Fourth  Picture  Coding  Confer¬ 
ence  to  be  held  on  22-24  January  1973.  Approximately  80  to  100 
research  workers  in  the  field  of  image  coding  are  expected  to  at¬ 
tend.  Topics  of  the  conference  include:  properties  of  the  human 
observer,  facsimile  coding,  intraframe  picture  coding,  interframe 
television  coding,  color  image  coding,  and  multispectral  image 
data  coding.  An  image  coding  contest  will  be  held  at  the  confer¬ 
ence  to  determine  the  best  image  coding  algorithms  and  systems 
for  monochrome  and  color  images. 
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3.  Image  Coding  Projects 

The  goal  of  the  image  processing  research  study  -  to  develop 
digital  image  coding  techniques  capable  of  representing  monochrome  and 
color  images  with  a  minimal  number  of  code  bits  -  is  not  likely  to  be 
satisfied,  at  least  in  the  near  term,  by  a  single  coding  system.  Accordingly 
the  research  effort  has  been  directed  toward  a  wide  variety  of  coding 
methods  designed  for  various  applications.  The  results  of  the  research 
study  during  the  past  six  months  are  summarized  here  and  presented  in 
greater  detail  in  the  subsequent  subsections. 

A  study  on  quantization  techniques  for  natural  color  images  is  the 
subject  of  the  first  report.  Using  a  quantitative  measure  of  color  error, 
previously  developed,  it  was  found  that  the  minimum  error  color  repre¬ 
sentation  using  standard  PCM  coding  techniques  is  the  red,  green,  blue 
representation.  Other  coordinate  systems  such  as  the  YIQ  system  employed 
for  commercial  television  transmission  lead  to  significantly  greater  quanti¬ 
zation  error  unless  a  complicated  quantizer  is  employed. 

The  study  on  deltamodulation  and  differential  PCM  coders  for  mono¬ 
chrome  and  color  images  has  led  to  worthwhile  results.  In  particular, 
simulation  of  a  dual  mode  coderwhich  uses  deltamodulation  for  low  detail 
image  regions  and  differential  PCM  for  high  detail  image  areas  has  indicated 
that  subjectively  i>ood  results  can  be  obtained  with  about  1.  5  to  Z.  0  bits /pixel. 
Such  a  coder  could  be  implemented  in  real  time  with  a  relatively  modest  system 

The  research  work  on  transform  coding  systems  has  continued  with 
the  development  of  the  slant  transform  coder.  The  slant  transform  possesses 
basis  functions  (expansion  waveforms)  that  more  closely  approximate  natural 
imagery  than  other  more  commonly  employed  transforms.  Analysis  and 
supporting  experimentation  has  indicated  that  the  Slant  transform  is  capable 
of  image  coding  with  about  1.0  to  Z.  0  bits/pixel.  Another  recent  develop¬ 
ment  in  the  application  of  transforms  to  image  coding  is  the  transform 
differential  coding  method.  In  this  system  a  unitary  transform  is  used 
to  obtain  a  low  detail  estimate  of  the  image.  The  pixel  differences  be¬ 
tween  the  original  and  estimate  are  then  coded.  With  this  system  it  is 
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possible  to  code  images  with  as  low  as  3.  5  bits/pixel  with  no  error  or  with 
about  1.  5  to  2.  0  bils /pixel  with  a  scarcely  noticeable  error.  It  now  appears 
that  these  transform  coding  systems  can  be  implemented  in  real  time  using 
acoustic  surface  wave  delay  line  implementation. 

A  problem  common  to  almost  all  types  of  image  coding  systems, 
is  the  efficient  coding  of  variables,  e.  g.  image  samples,  differential  PCM 
samples,  or  image  transform  coefficients.  A  comprehensive  comparison 
of  several  quantization  techniques  has  been  'ompleted. 

In  the  final  research  study  reported  upon  in  this  section,  some  pre¬ 
liminary  results  are  presented  on  a  new  coding  process  called  Universal 
Coding.  In  universal  coding  systems,  coding  is  accomplished  without 
complete  knowledge  of  the  source  statistics,  e.  g.  statistically  coding  an 
image  without  complete  knowledge  of  the  pixel  correlation.  Results  of 
this  theoretical  investigation  should  be  useful  in  establishing  performance 
bounds  on  image  coders  operating  with  incomplete  statistical  knowledge 
of  the  image  data. 

3.  1  Color  Image  Quantization 

Anil  K.  Jain 
■William  K.  Pratt 

When  color  images  are  converted  to  digital  form  for  computer  pro¬ 
cessing  or  digital  communication,  consideration  should  be  given  tc  the 
effects  of  quantization  errors.  Quantization  errors  can  cause  luminance, 
hue,  and  saturation  shifts  of  a  color  that  result  in  a  significant  image 
degradation.  In  quantization  of  color  images,  the  quantization  error 
depends  on:  (a)  the  relative  number  of  bits  assigned  to  the  different 
coordinates;  and  (b)  the  distribution  of  the  image  colors  in  the  color  space 
corresponding  to  the  respective  color  coordinate  system. 

Figure  1  contains  a  general  block  diagram  for  a  color  image  quanti¬ 
zation  system.  A  source  image  described  by  source  tristimulus  values 
R,  G,B  is  converted  to  three  components  Xj,X2,x^  which  are  then  quantized. 
The  quantized  components  x^.x^.x^,  are  then  converted  back  to  the  ori- 
ginal  color  coordinate  system  producing  the  tristimulus  values  R,  Q  B  . 
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Figure  3.1--1  Color  Image  Quantization  Model 
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The  quantizer  in  Figure  1  effectively  partitions  the  color  space  of  the 
color  coordinates  into  quantization  cells  and  assigns  a  single  color 

value  to  all  colors  with  a  cell.  To  be  most  efficient,  the  three  color  com¬ 
ponents  should  be  quantized  jointly.  However,  implementation 

problems  often  dicta.te  separate  quantization  of  the  color  components.  In 
such  a  system  individually  quantized  over  their  maximum 

ranges.  In  effect,  the  physical  color  solid  is  enclosed  in  a  rectangular  solid 
which  is  then  divided  into  rectangular  quantization  cells.  Many  of  these  quanti¬ 
zation  cells  will  lie  outside  the  physical  color  solid,  and  will  be  wasted. 

In  the  present  analysis  it  will  be  assumed  that  each  color  component  is 
linearly  quantized  over  its  maximum  range  into  2^^  levels  where  n.  represents 
the  number  of  bits  assigned  to  the  component  x..  The  total  number  oi  bits  al¬ 
lotted  for  the  coding  is  fixed  at 


N  =  n^  +  n^  + 


(1) 


Let  b.  represent  the  upper  bound  of  x.  and  a.  the  lower  bound.  Then  each 
1  11 

quantization  cell  has  dimension 


e 

1 


(2) 


Thus  the  coordinates  of  the  quantized  color  become 


y.  =  X.  ±  e.  (3) 

1  11 

subject  to  the  conditions 

a.  ^  y.  ^b.  ,  i  =  1,2,3  (4) 

111 

Observe  that  the  values  of  y  always  lie  within  the  smallest  rectangular 
solid  enclosing  the  color  solid  for  the  given  color  coordinate  system. 

The  total  error  due  to  quantization  can  be  specified  by  the  geodesic 
color  distance  [l"]  between  x.  and  y.  as  given  by 
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=  1, 


(yj 


(x.) 


1 
n  2 


E 

L 


C„  dx.dx, 
3^  J  ^ 


(5) 


where  the  coefficients  C  represent  the  sensitivity  of  average  human  per- 

th  th 

ception  of  differences  in  the  j  and  in  the  k  coordinates  of  a  color  coordinate 
system.  Using  equation  (5)  it  is  possible  to  evaluate  the  total  error 
due  to  a  given  bit  assignment.  Thus,  for  a  given  number  of  bits,  N,  it  is 
possible  to  find  the  optimal  bit  assignment  for  each  color.  However,  for  the 
coding  of  color  images,  one  would  like  to  determine  the  optimal  bit  assign¬ 
ments  for  all  the  colors  in  a  given  coordinate  system.  In  this  study  con¬ 
sideration  is  given  to  the  problem  of  determining  non-inferior  bit  assign¬ 
ments  for  the  N.  T.S.C.  receiver  primary  major  colors  in  the  following 
coordinate  systems:  i)  C.I.E.  Uniform  Chromaticity  scale;  UVW  and 
uvY ;  ii)  N.  T.S.C.  transmission  system,  YIQ;  iii)  N.  T.S.C.  receiver 
phosphor  system,  RGB  and  rgY.  Choosing  N  =  12,  equation  (5)  was  used 
(via  a  color  distance  routine  [l])  to  determine  the  maximum  value  of  S  ( over 
all  possible  y),  for  a  fixed  color  and  a  fixed  bit  assignment  n^.n^.n^  such 
that 


.  ^  6  ,  i  =  1,2,3 

(6) 

N  =  12 

(7) 

This  was  repeated  for  the  seven  major  colors  and  all  possible  bit  assign¬ 
ments  satisfying  (6)  and  (7). 

The  bit  assignments  giving  relatively  large  values  of  S  were  considered 
inferior  and  thus  ignored.  Table  I  for  example,  shows  the  non-inferior  set 
of  bit  assignments  for  various  colors  and  in  the  RGB  system  (see 
reference  [z]  for  more  results). 

Quantization  of  a  256  by  256  color  image  of  a  girl  [3],  was  done  with 
bit  assignments  (4,4,4),  (5,4,3)  and  (3,5,4)  in  the  RGB,  YIQ  and  UVW 
systems.  The  color  shifts  in  the  RGB  system  were  found  to  be  smallest 
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and  were  in  relative  agreement  with  the  error  values  of  Table  I.  In  the  UVW 
and  YIQ  systems,  relatively  large  amount  of  color  shifts  were  observed 
(compared  to  the  values  obtained  computationally  [z]  )  and  were  found  to  be 
very  sensitive  to  the  quantization  range  of  each  coordinate.  These  large 
color  shifts  occurred  because  the  actual  ranges  of  the  I,  Q  and  U,  W  com¬ 
ponents  of  the  image  were  much  smaller  than  the  corresponding  ranges  of 
these  coordinates  in  their  color  solids,  and  also  a  significant  nvimber  of 
quantization  cells  in  the  rectangular  solid  enclosing  the  actual  color  solid 
lie  outside  that  color  solid.  The  conclusion  from  this  study  is  that  for  a 
fixed  number  of  quantization  cells,  the  cells  should  be  packed  into  the 
volume  occupied  by  the  color  solid.  Such  a  quantization  procedure,  however, 
is  difficult  to  implement  except  for  the  RGB  coordinate  system. 


TABLE  I.  Quantization  errors  in  the  RGB  system  for  different 
non-inferior  bit  assignments. 


Bit 

Assignment 

Red 

Green 

'  ■  '  ■■1 

Blue 

Cyan 

Magenta 

Y  ellow 

White 

4,4,4. 

15.4 

7.0 

20.  6 

4.5 

8.  2 

5.0 

3.4 

4,5,3 

21.0 

10.3 

1.5 

5.  1 

3.6 

7.  8 

3.5 

5,3,4 

22.  8 

6.7 

2.9 

5.0 

13.  8 

0.8 

6.  0 

5,4,3 

23.2 

9.3 

13.9 

3.9 

7.  2 

7.0 

4.5 

3,5,4 

12.9 

11.5 

8.5 

0.  3 

7.7 

5.7 

3,4,5 

12.  3 

3.3 

8.  3 

0.  2 

6.8 

5.7 
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3.2  Delta  modulation  and  Differential  PCM  Coding  of  Monochrome  and 

Color  Signals 

A.  Habibi 

Deltamodulation  and  Differential  PCM  image  coding  systems  offer 
substantial  promise  for  real  time  television  coding  applications  because  of 
their  relatively  modest  implementation  requirements.  The  following  sec¬ 
tions  present  an  introductory  discussion  of  these  systems  and  descriptions 
of  their  application  to  the  coding  of  monochrome  and  color  images. 

DPCM  Systems.  A  differential  PCM  image  coder,  as  shown  in 
figure  Ij  consists  of  a  linear  predictor  that  estimates  the  incoming  pixel 
value,  a  subtractor  that  subtracts  the  estimated  pixel  value  from  its  actual 
value  to  form  a  differential  signal,  and  a  quantizer  that  maps  the  differen¬ 
tial  signal  into  a  set  of  discrete  levels.  An  efficient  predictor  for  video 
signals  is  one  that  exploits  the  correlations  of  the  signal  in  all  spatial 
directions.  Experimental  results  have  indicated  that  employing  a  third 
order  predictor  results  in  a  substantial  improvement  over  the  first  order 
predictor  [  1  ].  Further  increases  in  the  prediction  order  are  usually  not 
warranted. 

Deltamodulation  Systems.  In  a  deltamodulation  system  each  image 
sample  is  compared  to  its  estimate  and  a  positive  or  negative  signal  is 
produced  depending  on  the  comparative  amplitude  of  the  difference.  A 
block  diagram  of  a  delta  modulator  is  shown  in  figure  2.  The  output  of  the 
comparator  is  multiplied  by  a  constant  in  the  feedback  loop  and  is  used  as 
an  input  signal  to  an  integrator  whose  output  is  the  estimate  of  the  incoming 
signal.  The  size  of  the  constant  (step  size)  is  an  important  parameter  in 
the  encoder,  A  large  step  size  degrades  the  encoded  signal  when  the  signal 
is  changing  smoothly  by  causing  a  large  granular  noise;  a  small  step  size 
limits  the  ability  of  the  encoder  to  build  up  to  large  signal  change.  This  is 
avoided  by  introducing  some  memory  in  the  quantizer  that  considers  the 
polarity  of  the  previous  binary  digits  at  the  output  of  the  quantizer  and 
changes  the  step  size  accordingly.  A  small  step  size  for  a  slowly  varying 
signal  and  a  large  step  size  for  a  fast  varying  signal  enables  the  adaptive 
deltamodulator  to  track  the  input  signal  more  accurately.  Abate  [2]  has 
developed  a  workable  system  that  uses  step  sizes  of  1,  2,  and  4. 
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Figure  3.3-2  Delta  Modulation  Image  Coder 


Dual  Mode  Deltamodulation-DPCM  Systems*.  The  combination  of 
the  deltamodulation  and  the  DPCM  coding  techniques  has  been  found  to  pro¬ 
vide  an  improvement  in  information  transfer  rate  for  a  given  picture  quality 
for  monochrome  images  when  compared  to  either  system  separately.  Use 
of  DPCM  alone  is  inefficient  because  the  codeword  output  symbols  are  not 
typically  generated  with  equal  probability.  Use  of  deltamodulation  alone  is 
inefficient  because  of  quickly  changing  amplitude.  The  dual  mode  coder 
combines  the  best  of  both  modes  (deltamodulation  for  slowly  changing  regions 
and  DPCM  for  quickly  changing  regions)  and  potentially  achieves  a  bit  rate 
reduction  while  preserving  picture  quality.  The  simulation  results  con¬ 
firm  this  hypothesis. 

Figure  3  presents  a  pictorial  view  of  the  principal  somponents 
of  the  dual  mode  coder  system.  The  position  of  switches  SI  and  S2  determine, 
the  current  coding  mode.  The  switching  is  perform.ed  automatically  by  the 
system.  This  system  is  similar  to  one  considered  by  Frei  et  al  [sjwith  the 
difference  that  here  a  third  order  predictor  is  used  for  the  DPCM  system. 

The  quantizer  for  deltamodulation  is  uniform,  with  two  levels.  The 
deltamodulation  reconstruction  level  Q  determines  the  extent  of  mode  switch¬ 
ing  in  the  system.  The  DPCM  quantizer  is  nonuniform  with  8  levels. 

Prior  to  processing,  the  number  of  pixels  per  line  is  expanded  by  a 
factor  of  three.  The  new  pixels  are  interpolated  from  the  original  data. 

The  deltamodulation  coder  uses  the  expanded  lines  w'here  the  DPCM  coder 
uses  only  the  original  data.  During  the  deltamodulation  mode,  groups  of 
three  pixels  are  coded.  The  sequences  of  three  bits  produced  by  the  quan¬ 
tizer  are  then  compressed  to  one  bit  according  to  a  majority  decision.  If 
all  three  bits  are  identical  they  are  transmitted  uncompressed;  the  coder  also 
switches  to  DPCM.  During  the  DPCM  mode  each  of  the  original  pixels  is 
coded  for  transmission.  When  two  successive  transmissions  are  at  the 
lowest  quantization  level,  but  alternate  sign  (idling  condition),  the  coder 
switches  back  to  the  delta  mode  since  such  a  condition  indicates  the  resump¬ 
tion  of  a  smooth  picture  region.  One  drawback  of  this  dual  mode  coding  is 
the  generation  of  data  at  an  uneven  rate  which  requires  buffering  the  data 
prior  to  its  transmission. 

*  The  results  pertaining  to  the  combination  of  DPCM  and  deltamodulation  coder 
were  obtained  by  Thomas  Middleton  a  graduate  student  in  Electrical  Engineering 
Department  of  the  University  of  Sovithern  California. 
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Figure  3.Z-3  Dual  mode  deltamodulation  -  DPCM  image  coder 


Photographs  of  the  original  and  processed  pictures  appear  in 
figures  4  and  5,  Results  are  quite  encouraging  and  support  the  results 
reported  by  Frei  [3].  All  of  the  simulation  run  parameters  were 
trained  on  the  tank  picture.  With  the  delta  reconstruction  threshold  set 
at  Q=4,  0,  the  tank  showed  only  slight  degradation  coded  at  1, 72  bits/ 
pixel  (Huffman  DPCM  Coding).  Using  these  same  parameters,  the  girl 
was  coded  in  1.15  bits /pixel,  suffering  slightly  more  degradation.  With 
Q=5.  0  and  a  bit  rate  of  1.61  (3  bit  DPCM  coding)  the  tank  showed  slightly 
higher  contouring  but  still  a  good  picture.  At  Q=10.  0  and  1.16  bits /pixel, 
the  contouring  was  more  pronounced.  Thus  an  acceptable  picture  can 
probably  still  be  coded  at  about  1. 5  bits /pixel,  the  rate  achieved  by 
Frei  [3  ]  . 

Based  on  these  results  the  dual  mode  coder  appears  to  be  a 
workable  scheme,  particularly  in  a  noiseless  environment  such  as  stor¬ 
ing  video  data  digitally.  In  a  noisy  transmission  channel  coder  system 
it  is  not  as  attractive  because  redundancy  must  be  added  to  keep  per¬ 
formance  at  a  higher  level.  Differential  coding  schemes  are  unstable 
in  the  presence  of  noise  and  do  not  recover  until  reinitialization  such  as 
at  the  beginning  of  a  line. 

Comparison  of  DPCM  and  Pelt  Modulation  Systems  for  Color  Images. 
Both  the  encoder  and  the  decoder  of  a  DPCM  and  the  delta  modulation  sys¬ 
tems  were  simulated  on  a  digital  computer  to  code  the  luminance  and  the 
chromaticity  components  of  a  color  video  signal.  The  luminance  component 
of  this  video  signal  is  coded  by  the  DPCM  and  the  delta  modulation  systems 
at  various  bit  rates.  The  results  are  used  to  evaluate  the  performance  of 
each  coding  system  and  compare  their  performances. 

The  mean  square  values  of  the  coding  errors  (MSE)  are  plotted  in 
figure  6.  The  performance  of  both  the  first  and  the  third  order  DPCM 
systems  improve  at  about  the  same  rate  by  increasing  the  number  of  quanti¬ 
zation  levels.  The  first  order  DPCM  system  with  a  two-level  quantizer  is 
essentially  the  same  as  the  simple  delta  modulator  operating  at  one  bit  per 
picture  element.  The  adaptive  delta  mod\alator,  taking  advantage  of  the 
feedback  properties,  outperforms  both  of  these  systems  at  coding  rates  of 
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(b)  DPCM,  2.66 


10,  1.16  bits /pixel 


(i,  Dual 

Q  -  'I*  !•  18  bits/pixel 


Figure  3.3-4  Dual  mode  image  coding  using  "TANK”  picture 


f 


(a)  Original 


(b)  Dual 

Q  =  4,  1.16  bits/pixel 


(c)  lhaal 

Q  =  6,  1.07  bits/pixel 


Figure  3.3-5  Simulation  results  using  "GIRL"  picture 

-15- 


|uauJ9|3  sjnioy/s^jg  uj  040^  d|dq 


-16- 


Figure  3.2-6  Coding  Mean  Squared  Error  at  Various  Data  Rates  for  the  Normalized  Picture 


one  and  two  binary  digits  per  picture  element.  However,  the  first  order 
DPCM  system  achieves  a  better  coding  capability  at  higher  bit  rates.  As 
shown  in  figure  6  the  performance  of  the  third  order  BPCM  system  is 
better  than  the  performance  of  the  delta  modulators  at  all  bit  rates. 

The  results  of  the  delta  modulators  at  bit  rates  higher  than  one 
bit  per  picture  element  requires  sampling  the  analog  data  at  a  higher 
sampling  rate.  The  previously  sampled  data  was  interpolated  to  generate 
a  larger  number  of  samples.  The  linear  interpolation  used  to  double, 
triple,  and  quadruple  the  number  of  samples  introduces  higher  frequency 
errors  and  other  undesirable  distortions.  The  higher  frequency  error 
is  eliminated  in  practice  by  low  pass  filtering  the  coded  signal.  In  cal¬ 
culating  the  coding  error  for  the  simulated  systems  the  effect  of  the 
high  frequency  error  was  eliminated  by  two  different  techniques  (see  [4]). 
The  solid  lines  on  figure  6  corresponds  to  complete  elimination  of  these 
errors,  and  the  curves  shown  by  dashed  lines  correspond  to  partial 
elimination  of  the  high  frequency  errors. 

The  coded  pictures  corresponding  to  some  of  the  points  on  figure 
6  are  shown  on  figure  7.  The  degradation  in  the  picture  coded  by  the 
third  order  DPCM  system  at  a  bit  rate  of  2  bits  per  picture  element  is 
not  noticeable.  Using  the  other  scherues  a  rate  of  3  bits  per  picture 
element  is  needed  to  obtain  similar  results. 

The  color  video  signal  is  very  insensitive  to  degradations  in  the 
chromaticity  components,  thus  these  signals  can  be  coded  using  only  a 
fraction  of  binary  digits  per  signal  element.  Figure  8  is  a  display  of  the 
I  and  Q  signals  coded  by  the  adaptive  delta  modulator  using  1/4  bit  per 
pixel.  To  achieve  this  bit  rate  both  the  I  and  Q  signals  are  subsampled 
to  reduce  the  resolution  in  both  spatial  directions  by  a  factor  of  two. 

This  signal  is  then  coded  and  the  result  is  linearly  interpolated  to  in¬ 
crease  the  number  of  picture  elements  to  256  by  256.  Combining  these 
signals  with  the  coded  luminance  component  results  in  coded  video 
signal  which  does  not  show  any  noticeable  degradation.  The  total  bit 
rate  using  the  luminance  component  coded  by  the  third  order  DPCM 
system  is  2.  5  bits  per  picture  element.  The  bit  rate  corresponding 
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(a)  Original  (b)  3rd  Order  DPCM, 

2  bits/pixel 


(c)  Adaptive  Delta,  3  bits/pixel  (d)  Simple  Delta,  3  bits/pixel 


Figure  3.2-7  Lumination  of  a  Color  Signal  Coded  by  DPCM, 
Adaptive  Delta,  and  Simple  Delta  Modulators. 
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i 


(a)  Original  I 


(c)  I,  Adaptive  Delta 
0.  25  bits/pixel 


(b)  Original  Q 


■  - 


■  1 


(d)  Q,  Adaptive  Delta 
0.  25  bits/pixel 


-8.  Original  and  Coded  Chromaticity  Components  of  the  Color  Video  Signal 
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to  the  luminance  signal  coded  by  other  schemes  is  3.  5  bits  per  picture 
element. 
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3.3  Slant  Transform  Image  Coding  of  Monochrome  Images 
Weii-Hsiung  Chen  and  William  K.  Pratt 

In  transform  image  coding  using  separable  unitary  transforms 
the  two  dimensional  transform  F(u,  v)  of  an  image  f(x,  y)  may  be  written 
as 

[F(u,  v)]=  [a  ]  [f(x,  y)  ]  [a  ] 

where  [A]is  a  matrix  whose  rows  are  basis  functions  of  the  transforma¬ 
tion.  An  efficient  transformation  for  image  coding  is  one  for  which  the 
image  energy  is  redistributed  to  a  relatively  small  number  of  transform 
coefficients  F(u,v).  The  remaining  coefficients  can  then  be  quantized  with 
a  small  number  of  quantization  levels  without  significiant  error.  The  energy 
compaction  of  a  unitary  transform  will  be  "high"  if  the  basis  functions  of 
the  transform  "resemble"  typical  lines  of  an  image.  In  most  natural  images 
a  large  number  of  lines  or  line  segments  are  of  nearly  constant  brightness. 
Also  many  lines  and  line  segments  either  increase  or  decrease  in  brightness 
over  their  length  in  a  nearly  linear  fashion.  Thus,  one  of  the  basis  func¬ 
tions  should  be  a  constant  and  another  a  discrete  sawtooth  or  slant  wave¬ 
form.  With  this  thought  in  mind,  a  fam.ily  of  basis  functions  containing  a 
flat  and  a  sawtooth  waveform  member  have  been  developed.  The  remaining 
basis  functions  have  been  chosen  so  that  the  sequency  (number  of  zero  cross- 
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ings)  of  each  basis  function  is  equal  to  its  row  number  minus  one.  Also, 
by  design,  the  resulting  set  of  basis  functions,  called  the  Slant  trans¬ 
form,  possesses  a  fast  computation  algorithm.  The  algorithm  is  des¬ 
cribed  in  detail  in  ref.  [l  ].  Figure  1  is  a  superimposed  plot  of  the 
Slant  and  Hadamard  transform  basis  functions  for  a  vector  length  of  16. 

An  intensive  study  has  been  performed  on  techniques  for  efficiently 
quantizing  transform  domain  samples  for  the  commonly  employed  image 
transforms.  An  optimum  strategy  from  the  standpoint  of  minimizing  the 
mean  square  error  between  the  original  and  the  reconstructed  images  has 
been  found.  First,  the  number  of  bits  assigned  to  each  transform  coef¬ 
ficient  is  made  proportional  to  the  logarithm  of  the  variance  of  the  coef¬ 
ficient.  The  quantizer  levels  are  then  selected  according  to  the  Max 
quantization  rule.  Figure  2  contains  a  plot  of  the  mean  square  error  ob¬ 
tained  with  several  image  transforms  as  a  function  of  the  image  block 
size  processed  for  image  coding  with  an  average  of  1.  5  bits /pixel.  In 
this  example,  the  image  field  was  modelled  as  a  first  order  Markov  pro¬ 
cess  with  row  and  column  correlation  factors  of  0.  95.  From  this  figure 
it  is  seen  that  the  Slant  transform  is  much  superior  to  the  Hadamard 
transform  and  closely  approaches  the  optimal  performance  of  the  Karhunen- 
Loeve  transform.  Also,  the  curve  indicates  that  the  mean  square  error 
does  not  decrease  much  after  an  image  block  size  of  16  by  16  is  reached. 

Figure  3  contains  original  eight  bit,  256  by  256  pixel  images  and 
their  reconstructions  coded  with  an  average  of  1.  5  bits /pixel  using  the 
quantization  strategy  described  above.  Subjectively,  the  reconstructions 
appear  to  be  of  good  quality  with  little  visible  degradation. 

Further  research  on  slant  transform  image  coding  is  underway. 

The  effects  of  channel  errors  are  being  studied,  and  the  Slant  transform 
is  also  being  applied  to  color  images. 
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Figure  3.3-2  Mean  Square  Error  Performance  of  Transform  Image  Coders 
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3.4  Transform  Differential  Coding  of  Monochrome  Images 
Faramarz  Bavarian  and  William  K.  Pratt 

In  conventional  transform  image  coding  a  surprisingly  recog- 
nizeable  monochrome  image  can  be  reconstructed  using  about  5-10% 
of  the  transform  coefficients  with  the  largest  magnitude.  The  recon¬ 
structed  image  generally  possesses  the  correct  overall  grey  scale 
shading,  but  lacks  high  frequency  detail.  One  could  consider  this 
reconstructed  image  to  be  a  first  order  estimate  of  the  original  image 
in  the  t;ense  that  the  feedback  signal  in  a  differential  PCM  image  coder 
is  an  estimate  of  the  scanned  pixel  values.  This  idea  has  led  to  the 
development  of  a  coding  technique  called  transform  differential  image 
coding  as  shown  in  Figure  1.  In  this  system  an  image  block  undergoes 
a  two  dimensional  unitary  transform.  The  transform  coefficients  are 
grossly  quantized,  coded,  and  transmitted.  The  quantized  transform 
coefficients  are  then  subjected  to  an  inverse  two  dimensional  transfor¬ 
mation.  The  resultant  image  estimate  is  then  subtracted  pixel  by 
pixel  from  the  original  image.  The  pixel  difference  signals  are  quantized, 
coded,  and  transmitted.  At  the  receiver  the  pixel  difference  signals  are 
added  to  the  image  estimate  to  obtain  a  reconstruction  of  the  image  block. 

A  bandwidth  reduction  can  be  obtained  with  the  transform  differential  image 
coding  by  two  coding  techniques.  One  of  the  techniques  is  information  pre¬ 
serving  and  yields  a  moderate  compression  ratio,  and  the  other  technique 
permits  a  small  reconstruction  error  in  order  to  achieve  a  greater  com¬ 
pression  ratio. 

Information  Preserving  Coder.  In  the  information  preserving  coder 
version  of  the  transform  differential  coder,  the  pixel  difference  signals  are 
quantized  to  the  same  grey  scale  resolution  as  the  pixels  of  the  original 
image  and  then  coded  using  a  variable  length  Huffman  coder.  Typically, 
an  average  of  about  2.  0  to  2.  5  bits/pixel  is  allotted  to  the  pixel  difference 
signals.  It  should  be  noted  that  the  quantization  process  on  the  transform 
coefficients  does  not  create  any  reconstruction  error  since  the  pixel  dif¬ 
ference  signals  are  perfectly  coded,  but  the  transform  coefficient  quantization 
process  does  affect  the  coding  performance  of  the  overall  system.  Figure  2 


-25- 


IMAGE 


COEFFICIENTS 
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(b)  Decoding  system 


Figure  3.4-1.  Transform  dUfferential  image  coding  system 


illustrates  entropy  of  the  pixel  difference  signal  as  a  function  of  the 
number  of  transform  coefficients  transmitted  for  a  typical  picture 
using  a  16  by  16  pixel  block  slant  transform.  Figure  2  also  contains 
plots  of  the  number  of  bits  allotted  to  the  transform  coefficients  and 
the  total  number  of  code  bits  required  for  the  image  using  an  optimum 
coder  for  the  pixel  differences.  The  maximum  pixel  length  has  been 
limited  to  a  difference  of  five  quantization  levels  in  this  example.  If 
a  pixel  difference  greater  than  five  occurs  the  actual  pixel  difference 
is  coded.  The  simulations  indicate  that  for  an  optimized  system,  an 
image  can  be  coded  with  an  average  of  about  3.5  to  4.0  bits /pixel. 

This  result  is  somewhat  disappointing  since  a  coding  of  about  4.0 
bits/pixel  can  be  achieved  simply  by  coding  the  pixel  differences  along 
an  image  line. 

Non-Information  Preserving  Coder.  In  the  non-information  pre¬ 
serving  coder,  the  pixel  difference  signal  is  quantized  using  a  nonlinear 
quantizer  with  a  relatively  small  number  of  levels.  The  quantizer  design 
follows  the  same  general  rules  as  for  a  differential  PCM  image  coder. 

In  an  initial  set  of  simulation  experiments  shown  in  Figure  3,  the  slant 
transform  coefficients  were  coded  with  an  average  of  2.  0  bits /pixel 
and  the  pixel  difference  signal  was  quantized  to  only  three  levels  and 
coded  with  1.0  bits/pixel.  The  resulting  reconstruction  of  3.0  bits/pixel 
is  quite  good.  It  appears  possible  to  approach  a  coding  of  1.5  bits /pixel 
with  this  technique  using  an  adaptive  quantizer.  Further  investigations 
along  this  line  are  underway. 

3.  5  Quantization  Error  and  Entropy  of  a  Single  Random  Variate 
A.  Habibi 

Most  practical  coding  systems  that  are  used  to  convert  analog 
signals  to  binary  digits  are  composed  of  three  serial  processors.  The 
first  attempts  to  process  the  input  data  to  form  a  set  of  uncorrelated 
samples  with  a  continuous  amplitude;  the  second,  a  memoryless  quantizer, 
maps  the  continuous -valued  signal  to  a  set  of  a  given  number  of  discrete 
levels;  and  finally  the  third  processor  converts  this  data  to  a  set  of  binary 
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digits  by  a  fixed- length  or  a  variable  length  coder.  The  use  of  the 
memoryless  quantizer  is  justified  if  the  signal  at  the  input  of  the  quantizer  is 
uncorrelated.  Examples  of  decor  relating  processors  are  the  sampler 
and  linear  predictor  in  differential  PCM  and  the  Karhunen-Loeve  trans¬ 
formation  in  the  transform  coding  schemes  [  1,  2].  The  variable  coding 
scheme  in  the  output  of  the  quantizer  is  employed  to  generate  binary 
data  at  an  average  rate  equal  to  the  entropy  of  the  discrete  data.  This 
rate  is,  in  general,  lower  than  the  rate  of  data  generated  by  a  fixed- 
length  coder  which  is  equal  to  the  base-two  logarithm  of  the  number  of 
quantization  levels.  Naturally,  increasing  the  data  rate  will  reduce  the 
degradation  between  the  original  signal  and  one  that  is  reconstructed  from 
the  binary  digits  by  inverse  operations.  This  relationship  is  bounded  by 
Shannons  rate  distortion  function  which  is  the  absolute  bound  for  any  cod¬ 
ing  system  [3,4]. 

In  this  report  attention  is  focused  on  the  operation  of  the  memory¬ 
less  quantizer.  The  criteria  of  optimality  are  both  the  mean  square  error 
and  the  entropy  of  the  discrete  signal.  Thus  the  system  which  performs 
closest  to  the  rate  distortion  function  is  the  ideal  system.  Using  this 
criterion  the  performance  of  some  existing  quantizers  in  coding  uncor- 
relat'^d  signals  of  various  distributions  is  studied.  Of  specific  interest 
to  image  coding  are  the  Max  quantizer,  the  instaneous  companding, 
quantizer  and  the  uniform  quantizer  for  input  signals  possessing  Gaussian 
or  double- sided  exponential  probability  density  functions.  These  pro¬ 
bability  densities  are  of  interest  because  of  their  occurrance  in  DPCM  and 
transform  coding  schemes  for  pictorial  data. 


Review  of  Quantization  Schemes.  A  quantizer  maps  a  sample  x,  a 
continuous  variable  with  a  continuous  probability  density  P^(^)  into  one  of 

a  finite  set  [y^ . y^^  }  of  rational  numbers.  This  mapping  follows  the 

procedure: 

if  X.  <x<x.^j  then  x  is  mapped  to  y.(j  =  1,  .  .  . ,  N) 

J  J  J 

The  x.{j=2,  .  .  .  ,  N)  are  the  transition  levels,  with  x^  and  greatest 

lower  bound  and  least  upper  bound,  respectively,  to  the  input  signal,  and 
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the  y.(j  =  l.  2,  .  .  .  ,N)  are  the  quantization  values.  Given  the  total  number  N, 

of  discrete  levels  allowed,  a  quantizer  is  specified  by  the  sets  of  transition 

levels  X.  and  the  reconstruction  values  y.. 

J  J 

In  his  search  for  a  quantizer  with  a  minimum  distortion,  Max  Cs] 
minimized  the  expected  value  of  an  arbitrary  function  of  error  between 
the  input  and  the  reconstructed  signals,  and  obtained  a  set  of  sufficient 
conditions.  For  a  square  function  of  error  the  sufficient  conditions  are 


y.  =  2x.  -  y. 

J  J  J-1 


j=2,  ...,N 


X  .  , 


-y.)  p  (x)dx=0 


j  1,...,N 


Note  that  this  technique  minimizes  the  mean  square  error  without  con¬ 
sidering  the  entropy  of  the  quantized  numbers.  Equations  (1)  and  (2)  can 
only  be  solved  by  iterated  techniques.  Various  schemes  of  approximating 
the  above  quantizer  have  been  suggested  that  use  easier  techniques  of 
computing  the  transition  and  the  reconstruction  values  [6-8] 

In  a  uniform  quantizer,  as  shown  in  figure  1,  both  the  transition 
and  the  reconstruction  levels  are  spaced  at  equal  intervals,  q.  For  a 
given  number  of  quantization  levels  N,  the  maximum  quantization  levels 
is  X={N-1)  q/2.  The  parameter  q  can  be  chosen  to  minimize  the  mean 
square  error  [S],  or  to  minimize  the  entropy  of  the  quantized  signal  or 
a  combination  of  the  two  in  order  to  code  at  a  rate  closest  to  the  Shannon 
rate  distortion  function  [4].  The  Instantaneous  companding  quantizer  is  a 
third  type  of  quantir-er  that,  in  a  sense,  combines  attractive  features  of 
the  other  two  typos.  This  system  uses  a  uniform  quantizer  that  is  pr  - 
ceeded  by  companding  of  the  input  signal.  The  type  of  compression  de¬ 
pends  on  the  probability  density  function  of  the  signal.  The  quantizer  is 
followed  by  an  expander  that  compensates  for  the  effect  of  the  compression 
by  expanding  the  quantized  signal.  The  system  is  shown  in  figure  2  where 
f(. )  and  g(.)  are  the  pre-and  the  post-quantization  mappings.  Algazi  [6] 
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Figure  3.5-1^  Uniform  quantizer 


Figure  3.5-2.  Block  diagram  of  the  instantaneous  companding  quantizer 


n 


Figure'3.  3>3,  The  analog  model  of  the  instantaneous  companding  quantizer 
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modeled  this  quantizer  by  the  linear  system  of  figure  3.  Under  the  as¬ 
sumptions  that  the  quantization  noise  is  small,  additive  and  is  uncor¬ 
related  with  the  signal,  Algazzi  showed  that  the  system  will  minimize 
the  mean  square  error  if  mapping  f{. )  is  given  by 

^  1 

n  - 

f{x)  =  Kj  I  [P^{x)  ]  -^  dx  +  (3) 

and 

g  {.)  =  f'N.)  (4) 

The  same  results  have  been  obtained  by  Panter  and  Dite  using  a  dif¬ 
ferent  approach  [9,  10].  Mappings  f{.  )  and  g{.  )  are  obtained  from  (3) 
and  (4)  for  a  given  probability  density  P^{x)  .  Then  the  uniform  quanti¬ 
zation  parameter  (for  a  given  number  of  quantization  levels)  is  chosen 
to  minimize  the  mean  square  error  or  to  operate  at  the  closest  point 
to  the  rate  distortion  function. 

Comparison  of  Quantization  Schemes.  The  Max,uniform,  and 
the  instaneous  companding  quantizers  discussed  in  the  previous  sections 
have  been  designed  for  two-sided  exponential  and  Gaussian  probability 
density  functions.  These  probability  density  functions  have  been  chosen 
because  of  their  importance  in  image  coding.  Equations  (1)  and  (2) 
were  then  solved  iteratively  for  the  two-sided  exponential  and  Gaussian 
probability  density  functions  to  obtain  the  quantization  error  and  the 
entropy  of  the  quantized  signals  for  various  number  of  quantization  levels. 
The  number  of  binary  digits  required  to  achieve  this  coding  error  is  the 
base-two  logarithm  of  the  number  of  quantization  levels  for  the  system 
that  uses  fixed- length  code  words.  When  variable  length  coding  is  em¬ 
ployed  the  number  of  binary  digits  is  equal  to  the  entropy  of  the  quantized 
signal.  Figures  4  and  5  show  the  bit  rate  versus  the  mean  square  error 
for  the  exponential  and  Gaussian  probability  densities  for  fixed-length 
codes  while  figures  6  and  7  show  a  similar  relationship  for  the  variable- 
length  coder.  In  figures  4  and  5  the  parameter  q  of  the  uniform  quantizer 
is  chosen  to  minimize  the  mean  square  error.  The  parameter  q  varies 
over  a  range  of  values  for  different  numbers  of  quantization  levels  in 
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it  rate  versus  the  mean  square  error  for  fixed-length  code  words- exponent;’ al  density 


Figure  3.  5-6.  Bit  rate  versus  the  mean  square  error  for  variable-length  coder -exponential  density 
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Figure  3.  5-7.  Bit  rate  versus  the  mean  square  error  for  variable -length  coder-Gaussian  density 


figures  6  and  7.  For  the  Gaussian  probability  density  function  the  tan¬ 
gent  to  these  curves  is  within  1/4  of  a  binary  digit  from  the  rate  dis¬ 
tortion  function,  and  for  the  two-sided  probability  density  it  is  within 
1/8  of  a  binary  digit.  Figures  4  and  5  indicate  that  for  a  fixed- length 
coder  the  efficiency  of  the  instantaneous  companding  quantizer  is  worse 
than  the  Meix  quantizer,  but  it  is  better  than  the  uniform  quantizer. 

Figures  6  and  7  show  that  the  instaneous  companding  quantizer  does  not 
perform  as  well  as  the  uniform  quantizer  using  a  variable  coding 
scheme,  but  ic  outperforms  the  Max  quantizer. 

From  the  results  above  it  is  concluded  that:  (1)  If  the  quantizer 
is  followed  by  a  variable  length  coding  scheme  (such  as  the  Huffman  Coder) 
then  the  best  results  will  be  obtained  by  using  a  uniform  quantizer.  The 
performance  is  not  very  sensitive  to  the  least  upper  bound  of  the  transition 
levels  X.  Good  coding  results  are  obtained  with  a  value  of  X  about  three 
times  of  the  standard  deviation  of  the  signal.  The  performance  of  the 
instantaneous  companding  quantizer  is  close  to  the  uniform  quantizer,  but 
it  remains  inferior  to  it  with  the  Max  quantizer  doing  worse  than  both 
schemes.  (2)  Using  fixed-length  code  words  the  performance  of  the 
Max  quantizer  is  superior  closely  followed  by  the  instantaneous  companding 
and  the  uniform  quantizers. 
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3.  6  Universal  Coding 
Lee  D.  Davisson 

Suppose  one  is  given  a  time  and  state  discrete  source  (e.g.,  a 
sampled  and  quantized  image)  to  block  encode  without  error  to  minimize 
the  average  code  word  length,  the  code  depending  on  some  unknown  para¬ 
meters  of  the  source  message  probabilities.  The  unknown  parameter 
couid  be  something  as  simple  as,  say,  the  probability  of  a  one  in  a  binary 
independent  letter  source  or  the  transition  probability  matrix  for  a  Markov 
model  of  an  image,  or  as  general  as  the  complete  set  of  message  proba¬ 
bilities  for  the  source.  Such  sources  are  called  composite. 

In  some  cases,  although  not  usually,  the  unknown  parameter  may 
itself  have  a  known  prior  distribution  so  that  all  ensemble  message  pro¬ 
babilities  are  known.  There  is  no  reason  to  suppose,  however,  that  cod¬ 
ing  for  the  whole  ensemble  will  be  efficient  for  the  actual  parameter  in  effect. 
One  of  the  important  results  of  this  study  is  to  establish  that  coding  for  the 
whole  ensemble  does,  under  certain  circumstances  have  a  very  important 
meaning  in  an  asymptotic  sense.  For  example,  consider  a  composite 
binary  source  whereby  the  probability,  6 ,  of  a  "one"  is  chosen  by  "nature" 
randomly  according  to  the  uniform  probability  law  on  [O,  1  ]and  then  fixed 
for  all  time.  The  user  is  not  told  the  8,  however.  Subsequently,  the 
source  produces  independent  letter  outputs  with  the  chosen,  but  unknown, 
probability.  The  probability,  given  9,  that  any  message  block  of  length  N 
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contains  n  ones  in  any  given  pattern  is  then 

(1) 

If  6 were  known,  a  per  letter  average  code  word  length  of 

H  (9)  =  -eiog^  9-  (1  -  8)  log^  (1  -  9)  (2) 

bits  would  be  required  to  block  encode  the  source.  Although  9 is  unknown, 
the  composite  source  output  probabilities  are  known.  In  fact,  the  pro¬ 
bability  of  any  given  message  containing  ones  in  N  outputs  in  simply 

y*  9"  (l-9)^~"d9  ^  1  (3) 

(Note  that  the  composite  source  does  not  produce  independent  letters. 

This  is  due  to  the  fact  that  past  observations  provide  information  about 
the  value  of  9  in  effect.  ) 

There  is  no  reason  to  suppose  that  the  use  of  the  probabilities  of  (3) 

to  design  a  variable  length  code  will  do  well  for  any  9  that  comes  along  in 

the  sense  of  approaching  (2)  by  using  large  blocks.  One  of  the  surprising 

results  of  this  study  is  that  the  probabilities  of  (3)  provide  codes 

as  good  asymptotically  as  the  probabilities  of  (1)  with  9  known  exactly. 

Consider  the  following  encoding  procedure.  It  is  seen  from  (1) 

that  all  sequences  with  the  same  number  of  ones  are  equally  probable. 

Send  the  message  in  two  parts.  Part  one  consists  of  a  fixed  length  code 

for  the  number  of  ones  using  at  most  log  (N+l)  +  l  bits.  Part  two  is  a  code 

/N\ 

indicating  the  location  of  the  ones  within  the  block,  one  of  (  )  possibilities. 

/N  \  ^ 

At  most  the  log  I  )  +1  bits,  a  variable  length  depending  on  the  number  of 

'  n  ' 

ones,  is  required  for  this  information.  The  per  letter  coding  rate  is  then 

i  log(N  +  l)+2)+log  (^)  ]  (4) 

The  first  part  goes  to  zero  as  N  The  expected  value  of  the  second 

part  must  converge  to  the  per  letter  entropy  of  equation  (2)  because  the  cod¬ 
ing  rate  cannot  be  below  the  source  entropy,  and,  at  the  same  time,  the 
second  part  of  the  code  is  optimal,  all  messages  conditioned  on  the  first 
part  being  equally  likely,  independent  of  9.  This  will  be  demonstrated  quan- 
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titatively  subsequently.  Note  that  a  comparison  of  (4)  and  (3)  shows  that 
coding  with  respect  to  the  whole  ensemble  with  6  uniform  on  [O,  1  ]  is 
asymptotically  optimum  for  this  example. 

Formulation  of  the  Universal  Coding  Problem.  Several  measures 
of  optimal  universal  coding  effectiveness  can  be  proposed  in  terms  of  a 
source  block  output  random  variable,  an  unknown  parameter  _0  . 

As  they  all  involve  differences  between  the  average  code  word  length, 
denoted  i  (3^^  |  |_)  and  the  entropy  per  block  symbol,  H(  |  9),  the 
measure  is  a  form  of  redundancy.  In  each  case  the  infimum  is  taken 
over  the  class  of  uniquely  decipherable  codes  on  block  N.  They  are 
listed  in  order  of  ascending  redundancy: 

Definition  1.  Weighted  Redundancy 


w 


lim  inf 
N  C 


/ 


d»(8) 


N 


N 


where  w  (^)  is  a  probability  measure. 

for  which  R  approaches  zero 
w 

Definition  2.  Maximin  Redundancy 


A  code  for  which  R  approaches  zero  is  called  weighted  universal. 


^  =  lim  sup  inf 
N  -»«  weW  C 


N 


f  £{X 

Ja 


N 


■dw  (6 ) 


N 

If  ^  is  zero,  the  code  sequence  which  approaches  it  is  called  maximin- 
universal. 

Definition  3.  Minimax  Redundancy 


^  =  lim 

N  -♦< 


inf  sup 
a,  6  e  A 


N 


N  - 

If  the  minj  ax  redundancy  R  is  zero,  then  the  user  can  transmit  the  source 
output  at  a  te  per  symbol  arbitrarily  close  to  the  source  entropy  rate  uni¬ 
formly  over  all  values  ^  e  A  (unlike  the  other  forms  of  universality)  by  taking 
N  large  enough,  =  0  is  the  strongest  and  most  desirable  condition.  However, 
in  some  cases  it  is  too  much  to  hope  for.  If  a  sequence  of  codes  for  which 
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^  =  0  does  exist,  it  is  called  minimax- universal. 
It  can  be  shown  that: 


Theorem  1.  The  necessary  and  sufficient  condition  for  the  exis¬ 
tence  of  weighted  universal  codes  is  that  the  per  letter  average  mutual 
information  between  the  message  ensemble  {x.  }and  ^  be  zero. 

Theorem  2.  Meiximin  universal  codes  exist  if  and  only  if  the 
per  letter  channel  capacity  between  ^  and  {x^}  is  zero. 

Theorem  3.  A  necessary  and  sufficient  condition  for  the  exis¬ 
tence  of  minimax  universal  codes  is  that  there  exist  a  mixture  probability 
mass  function  p(Xj^)  for  which  the  conditional  mutual  information  be  zero 
uniformly  in  ^  and  N,  i.  e. 

lim  ^  sup  I(x^,  9^)  =  0 
N  0  eA 


The  determination  of  the  existence  or  non-existence  minimax  uni¬ 
versal  codes  is  mot  easily  made  by  searching  the  class  of  weighted-uni¬ 
versal  codes.  This  is  analogous  to  decision  theory  where  Bayes  rules  are 
more  easily  formed  than  minimax  rules  so  that  minimax  rules  are  con¬ 
structed  by  searching  the  class  of  admissible  Bayes  rules  by  finding  the 
Bayes  rule  with  respect  to  a  "best  guess"  least  fcivorable  prior. 

For  the  binary  example  presented  here: 


N 

5  i)  =  5  Z)  (H)  «"  (i-e)^'"iog[(N+i)(^) 

n=0 

=  log(N  +  l)  -  J_  .  ("entropy  of  a  binomial  random  j 
N  N  Lvariable  with  parameter  9 

<  log(N  f 1) 

N 

-  0 
N  -»® 

Hence,  the  code  constructed  by  taking  ©uniform  results  in  a  minirnax  uni¬ 
versal  code,  as  well  as  a  maximum  and  weighted  universal  code. 

Summary.  The  preceeding  has  charted  the  basic  concepts  of  uni¬ 
versal  coding  for  arbitrary  information  sources.  Studies  are  underway  to 
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apply  universal  coding  to  theoretical  image  models  and  to  digitized  images. 


3.7  Analog  Real  Time  Implementation  of  Transform  Image  Coders* 

Lloyd  R.  Welch 

The  invention  of  surface  wave  acoustic  devices  has  made  possible 
the  design  of  tapped  dely  line  filters  of  long  lengths.  (Filters  with  several 
hundred  delays  are  practical. )  Since  an  N  point  transform  can  be  thought 
of  as  N,  tapped  delay  line  filters,  the  thought  naturally  arises  as  to  the 
possibility  of  the  real  time  mechanization  of  transforms  using  surface 
wave  devices . 

Since  N  distinct  surface  wave  substrates  with  associate  driving 
circuits  would  be  expensive,  various  investigators  have  considered  mul¬ 
tiple  use  of  a  delay  line  by  connecting  the  taps  through  semiconductor 
switches.  A  switch  arrangement  is  shown  in  figure  1  when  the  transform 
coefficients  are  ^1,0.  The  number  of  such  taps  is  2N-1.  To  produce 


the  inner  product  of  the  signal  with  the  first  row  of  the  transform  matrix  at 

time  1,  through  are  appropriately  activated  while  for  k  >M  are 

_  ±  ±2 

all  off.  To  produce  the  inner  product  with  the  second  row  at  time  2,  S  to 


are  appropriately  activated  while  and  for  k  >N+1  are  all  off. 


It  is  somewhat  doubtful,  however,  that  the  signal  amplitude  out  of  the 


switches  can  be  controlled  well  enough  to  produce  an  accurate  summation. 
Furthermore,  when  the  transform  has  coefficients  with  amplitudes  other  than 
unity,  the  switches  cannot  be  controlled  well  enough  to  produce  these  am¬ 
plitudes  accurately. 

Recent  studies  by  material  science  investigators  have  shown  that 
it  is  not  difficult  to  lay  several  paths  on  a  substrate  driven  by  a  single 
transducer.  This  offers  an  alternative  configuration  where  each  row  of 
the  transform  matrix  corresponds  to  a  delay  line  channel  on  the  substrate 
with  its  own  tap  configuration.  Since  the  tap  apertures  can  be  precisely 
cut,  summation  coefficients  can  be  controlled  easily  to  less  than  one  percent. 
Several  problems  remain,  however.  Multichannel  devices  result  in  parallel 
output.  If  several  outputs  are  desired,  the  tap  configuration  on  consecutive 
delay  lines  must  be  offset  by  one  sample  time,  thus  the  total  amount  of  delay 

*  supported  in  part  by  the  Air  Force  Office  of  Scientific  Research 
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needed  on  the  substrate  is  still  2N-1  units  of  sample  time.  A  second  pro¬ 
blem  not  so  easily  solved  is  as  follows:  The  signal  on  the  substrate  is  g(t) 
cos  (2  TTft)  where  f  is  the  carrier  frequency  and  g(t)  is  either  the  signal 
to  be  processed  or  the  signal  with  a  bias  term  added.  In  either  case,  the 
tap  configuration  on  many  of  the  channels  will  be  orthogonal  to  the  bias 
term  and  the  output  g  cos  (  2nft)  may  have  a  negative  value  for  g  .  There¬ 
fore  a  simple  amplitude  detector  will  preserve  the  absolute  value  but  lose 
the  algebraic  sign.  The  method  used  to  overcome  this  problem  will  de¬ 
pend  on  what  further  analogue  digital  processing  will  be  done.  This  pro¬ 
blem  is  currently  under  investigation. 


4.  Image  Enhancement  and  Restoration  Projects 


Image  enhancement,  as  defined  here,  is  the  manipulation  of  pictoral 
data  for  improvement  in  presentation  to  a  human  viewer,  or  for  simplification 
in  format  for  subsequent  machine  processing.  In  image  enhancement  pro¬ 
cesses  there  is  no  overt  attempt  to  preserve  or  improve  fidelity.  Image 
restoration  processes  attempt  to  reconstruct  an  image  according  to  some 
pre -conceived  notion  of  what  the  image  should  have  been  if  there  had  been 
no  distortion  in  its  formation.  The  latter  problem  of  image  restoration  is 
amenable  to  hard  analysis,  while  the  former  problem  of  image  enhancement 
is  usually  subject  to  heuristic  approaches.  Image  enhancement  and  restoration 
research  topics  considered  during  the  past  six  months  are  summarized  below. 

The  first  three  research  tasks  deal  with  the  application  of  generalized 
Wiener  filtering  techniques  to  the  restoration  of  images  that  are  statistically 
described  in  terms  of  covariance  functions  of  the  original  image  and  the  noise 
or  interference.  In  the  first  report  consideration  is  given  to  restoration  of 
images  subjected  to  a  linear,  but  not  necessarily  space  invariant,  blurring 
and  additive  noise.  A  process  has  been  developed  for  sequential  row  and 
column  filtering  of  the  image  which  yields  substantial  improvements  in 
image  quality.  The  second  report  is  concerned  with  the  computat: onal 
aspects  of  generalized  Wiener  filtering  for  images  subjected  to  additive 
noise.  A  procedure  for  restricting  the  Wiener  filter  to  be  causal,  i.  e.  , 
only  operating  on  past  data  samples,  is  given.  It  is  demonstrated  that 
the  increase  in  reconstruction  error  under  a  causality  condition  is  small 
if  Hadamard  transform  preprocessing  of  the  data  is  employed.  The  third 
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report  deals  with  the  combination  of  the  scalar  Wiener  filtering  and  the 
homomorphic  filtering  concepts.  A  logarithmic  operation  is  performed 
on  the  image  to  separate  multiplicative  interference  from  an  image.  Linear 
Wiener  filtering  is  then  employed  to  minimize  the  interference,  followed  by 
an  exponentiation  operation  to  restore  the  image.  Significant  improvements 
in  the  detection  of  image  detail  are  demonstrated. 

The  following  two  reports  discuss  the  problem  of  image  restoration 
for  images  subjected  to  linear  space  variant  degradation.  It  is  shown  that, 
in  a  wide  class  of  problems,  the  space  variant  operator  can  be  cast  into  the 
cascade  of  an  invertible  geometric  operator,  a  linear  space  invariant  system, 
and  another  invertible  geometric  operator.  Such  systems  can  be  "undone" 
and  restored  by  conventional  space  invariant  techniques  such  as  scalar 
Wiener  filtering.  An  example  of  space  variant  restoration  is  given  for  an 
image  subjected  to  motion  during  film  exposure. 

The  next  report  outlines  some  beginning  attempts  at  image  restoration 
by  nonlinear  filtering.  It  is  pointed  out  that  conventional  linear  filtering 
often  results  in  non-physically  realizable  reconstructions,  i.  e.,  the 
representation  of  negative  intensities.  A  nonlinear  recursive  technique 
is  developed  for  image  restoration,  and  examples  of  its  performance  are 
presented. 

The  last  report  is  in  the  form  of  an  inserted  reprint  from  the  July  1972 
issue  of  the  IEEE  Spectrum.  This  paper  outlines  several  image  enhancement 
techniques  and  discusses  the  uses  of  pseudocolor  for  image  enhancement. 
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4.  1  Image  Deblurring  by  Generalized  Wiener  Filtering 
William  K.  Pratt 

The  concept  of  generalized  Wiener  filtering  previously  developed  [l^ 
for  the  restoration  of  images  degraded  by  additive  sensor  noise  can  be 
extended  to  compensate  for  spatial  blurring  of  the  image  as  well  as  noise. 


Image  Degradation  Model.  Figure  1  contains  a  block  diagram  of  a 
model  for  image  degradation.  The  ideal  continuous  image  is  subjected  to  a 
spatial  blurring  mechanism  which  is  assumed  linear,  but  not  necessarily 
spatially  invariant.  Next,  the  blurred  image  is  combined  with  additive 
noise  that  is  uncorrelated  with  the  image,  and  the  blurred  and  noisv  image  is 


then  spatially  sampled  and  quantized  for  digital  processing.  The  spatially 
sampled  image  function  can  be  expressed  as 


=  [fg(x,  y)  +  n{x,  y)]  s{x,  y) 


where  the  blurred  image  is  given  by  the  superposition  integral 


fg(x,y)  =  J  fj(a,  p)b{x,  y;a,  P)dadP 


with  b(x,y;a,P)  representing  the  impulse  response  of  the  physical  blur 
mechanism.  The  sampling  array 

M-1  N-1 

s(x,  y)  =  AxAy  E  E  j(x-mAx,  y-nAy)  (3) 

m=-M  n=-N 

is  assumed  to  be  composed  of  (M-N)  identical  sample  pulses  j(x,  y), 
arranged  in  a  grid  of  spacing  (Ax,  Ay).  The  image  at  the  sample  points 
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Figure  4.1-1  Model  for  image  degradation 
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Figure  4.  1-2.  Geometric  arrangement  of  image  sample  points 


(j  =  mAx,  k  =  nAy)  can  be  approximated  as  the  double  summation 


j+c  k-d 

fg{j,k)  =  fj{p,q)b{j,k;p,q)  +  n{j,k)  (4) 

p=j-c  q=k-d 

where  the  constants  (c,  d)  define  an  enclosing  boundary  for  the  blur 
mechanism  impulse  response.  Figure  2  illustrates  the  geometric  arrange¬ 
ment  of  sample  points  in  the  ideal  image,  blur  impulse  response,  and  the 
sampled  blurry  image.  It  should  be  noted  that  the  sample  points  outside 
the  physical  sampling  array  contribute  to  samples  of  the  blurred  image. 

In  general  the  blur  impulse  response  is  a  four  dimensional  function 
that  is  spatially  convolved  with  the  ideal  image.  If  the  blur  function  is 
separable  in  orthogonal  coordinates,  i.  e. 

b(j,k;p,q)  =  b^(j;p)bj^(k;q)  (5) 

then  the  rows  and  columns  of  the  image  array  can  be  processed  sequentially. 
In  this  case  eq.  (4)  can  be  expressed  in  matrix  form  as 

[fg(j,k)'l  =  [bj(j:p)l  [fj(p,q)]  [b^(k:q)]'^  +  [n(j,k)l  (6) 

The  development  of  the  generalized  Wiener  filter  for  image  deblurring 
and  noise  smoothing  presented  here  is  limited  to  the  sequential  processing 
of  rows  and  columns  of  an  image.  The  more  general  case  of  two  dimensional 
filtering  is  under  investigation. 

Sequential  Row  ?.nd  Column  Wiener  Filtering.  For  notational  simplicity 
let 
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^  =  P  X  1  column  vector  of  ideal  image  row  or  column 

^  =  M  X  1  column  vector  of  observed  image  row  or  column 

n_  =  M  X  1  column  vector  of  noise  along  image  row  or  column 

b  =  M  X  P  matrix  of  blur  impulse  response  along  image  row  or  column 

Then,  the  observed  image  vector  can  be  expressed  as 

^  =  b  ^  +  n  (7) 

In  the  generalized  Wiener  filtering  solution  to  the  image  restoration 
problem,  as  illustrated  in  Figure  3,  an  estimate  ^  of  the  ideal  image 
line  £  is  obtained  by  performing  a  unitary  transform  on  the  observed  image 
line  multiplying  the  transformed  data  vector  by  a  filter  matrix  G  ,  and 
performing  an  inverse  transform.  The  resulting  output  vector  is  then 

1  =  A'^GaAI  (8) 

The  reason  for  introducing  the  unitary  transform  A  is  to  decorrelate  the 
observed  data  vector  and  permit  subsequent  computational  simplifications 
in  the  filter  function  multiplication.  Transforms  that  have  been  considered 
include  the  Fourier,  Hadamard,  and  Karhunen- Loeve  transforms.  This 
subject  has  been  covered  in  Ref.  [l^  and  is  extended  in  the  following  section 
for  the  case  of  additive  sensor  noise,  but  no  image  blur.  In  developing  the 
optimum  filter  matrix  G  .  for  a  particular  unitary  transform  matrix  A,  it 
is  possible  to  compute  the  optimum  filter  matrix  for  the  identity  transform, 
and  then  perform  the  two  dimensional  transformation 
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(9) 


=  AG  A 
—A - 1 — 


The  optimum  Wiener  filter  matrix  is,  by  definition,  the  filter  matrix  that 
minimizes  the  mean  square  error  between  £  and  its  estimate  £.  However, 
the  estimate  will  only  be  used  over  the  vector  length  of  the  observed  data. 
Thus,  the  mean  square  error  may  be  written  as  the  trace  of  the  covariance 
matrix  of  the  error  over  the  data  vector  length. 


where  k  is  an  identity  matrix  of  dimension  M  embedded  in  the  center  of  a 
matrix  of  null  elements  of  dimension  P 


For  the  case  in  which  the  ideal  image  and  noise  are  uncorrelated,  the 
optimum  filter  matrix  is  found  to  be 
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(12) 


G  =  kC  b^fbC  +  C  1 

—I - S— - 8—  - 


-1 


and  the  resulting  mean  square  error  is 


(e)  .  =  kC  k  -  Gb  C  k 

min - 8—  — - s 


(13) 


In  the  limiting  case  in  which  there  is  no  blur,  the  blur  matrix  becomes 

f - P  - > 


/>■ 

M 

0 

_ 1 

^P/  2  — - M - )|f-P  /  2  • 


The  optimum  Wiener  filter 

G,  =  kC  k[k  C  k  +  C  1  (14) 

—I - s— - s—  — n 

then  reduces  to  the  form  derived  in  Ref.  [l]  for  no  blur.  If  there  is 
no  image  noise  the  optimum  Wiener  filter  becomes 


G^  =  k  C  b^[b  C  b  1 
—I - s— - s— 


Tn-l 


(15) 


Image  Restoration  Example.  Consider  the  restoration  process  when 
the  image  lines  are  statistically  described  by  a  first  order  process  with  a 
covariance  function 


C 


(16) 
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and  the  noise  is  white  with  a  covariance  function 

where  SNR  is  the  image  signal-to-noise  ratio.  The  blur  function  is 
assumed  to  be  Gaussian  in  form 

b{i,j}  =  exp{-(i-j)^/Bj^}  (18) 

Figure  4  is  a  plot  of  the  variance  of  each  pixel  of  the  image  line  estimate 
as  a  function  of  signal-to-noise  ratio  and  blur  for  an  image  line  section 
of  32  pixels  and  an  observation  of  16  pixels.  The  image  correlation 
factor  is  P  =  0.  8.  The  reason  for  utilizing  only  the  center  16  pixels  of 
the  32  pixel  estimate  is  readily  apparent  from  Figure  4. 

A  series  of  simulations  have  been  performed  to  evaluate  the  filtering 
process.  An  original  image  of  256  by  256  pixels  and  6  grey  levels 
wa&  filtered  over  its  entire  size  with  the  filter  function  of  eq.  (18)  using 
Fourier  domain  technique.  Numerically  generated,  zero  mean  Gaussian 
noise  was  ?,dded  to  each  pixel  value.  The  resultant  blurry  and  noisy  image 
was  then  requantized  to  64  grey  levels  to  serve  as  an  input  to  the  restor¬ 
ation  filter.  Figures  5a,  5c,  and  5e  are  photographs  of  blurred  and 
noisy  originals.  The  corresponding  restorations  are  shown  in  Figures  5b,  5d, 
and  5f.  It  is  clear  that  the  restoration  filter  has  resulted  in  a  substantial 
visual  improvement.  Also,  there  is  no  apparent  checkering  effect  in  the 
images  even  though  the  filtering  is  limited  to  16  X  16  pixel  blocks. 
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Figure  4.  1-4.  Variance  of  pixel  eatimatea 
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Figure  4.  1-5.  Generalized  Wiener  filtering  image  deblurring  examples 
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4.  2  Causal  Generalized  Wiener  Filtering  for  Image  Restoration 
Nelson  D.  A.  Mascarenhas 

Generalized  Wiener  filtering  computational  techniques  have  been 
used  with  excellen  results  in  one  and,  more  recently,  two  dimensional 
processing.  A  typical  example  of  the  last  case  is  the  problem  of  the 
restoration  of  an  image  corrupted  by  additive  noise  Cl].  It  has  also  been 
shown  that  suboptirnal  filtering  can  be  done,  with  reduced  computational 
demands,  without  severe  degradation  in  performance.  This  section  pre¬ 
sents  the  derivation  of  a  suboptirnal  Wiener  filter  satisfying  the  additional 
constraint  that  the  matrix  multiplication  that  performs  the  filtering  must 
be  causal  or  with  a  preassigned  delay.  Major  advantages  result  from  this 
restriction.  First,  the  number  of  operations  on  the  filter  matrix  multi¬ 
plication  can  be  considerably  decreased.  Second,  the  delay  in  starting 
this  operation  is  reduced.  This  also  leads  to  a  further  increase  on  the 
speed  of  the  overall  filtering  operation. 

Derivation  of  the  Suboptirnal  Filter.  The  discrete  causal  filter  for 
a  mean  square  error  criterion  was  first  derived  by  Friedland  [2].  The 
concept  is  extended  here  by  permitting  unitary  preprocessing  on  the  data. 
Figure  3  of  Section  4.  1  shows  the  block  diagram  of  a  generalized  Wiener 
filter  in  a  system  for  sequential  image  row  and  column  processing.  In 
this  system  £  is  the  input  (Mxl)  data  vector,  A  is  a  unitary  (MxM)  matrix 
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and  G  is  the  (MXM)  filter  matrix.  No  assumption  is  made  on  the 

stationarity  of  the  input  process,  but  it  is  assxuned,  for  simplicity,  that 

the  signal  and  noise  are  uncorrelated,  although  the  analysis  can  easily 

be  extended  to  the  other  case.  The  covariance  matrices  of  the  signal  ) 

and  noise  (C  )  are  assumed  known. 

— n 

The  design  of  the  filter  G .  is  based  on  minimizing  the  mean  square 


error 


e  =  Tr{  (s-s)(s*-s*)  } 


However,  since  A  is  a  unitary  transformation,  it  is  equivalent  to  minimize 


e  =  Tr{(S-S)(S*-S*)^} 


Observe  that  if  ^  is  a  stationary  process,  the  process  ^  may  be  nonstationary. 
Also,  if  the  matrix  A  represents  a  time  invariant  filter  (i.  e.  A(i,j)  =  A(i-j)) 
then  ^  will  be  stationary. 

The  optimum  noncausal  filter  has  been  derived  as  [l] 


Ga  = 

where 

- 1  .-T 

C  =  A  C  A  =  A  C  A- 
— S - s —  - s — 


C_  =  A  C  A"^  =  AC  A-"^ 
— - n —  - n — 


In  order  to  obtain  an  expression  for  G  in  the  suboptimal  filter,  first  consider 
the  case  where  no  del^  is  allowed  on  the  product  of  G^  by  This  restricts 


G .  to  be  a  lower  triangular  matrix.  A  direct  approach  to  obtain  G  in 
this  case  would  consist  in  imposing  the  orthogonality  between  the  error 

A 

(S-S)  and  the  data  (^)  vectors.  In  the  general,  non -stationary  case,  a 
set  of  1  +  2  +  •  •  •  +  M  =  M(M+l)/2  linear  equations  would  result  in  the 
same  number  of  elements  of  G  below  or  on  the  main  diagonal.  In  the 
stationary  case  this  number  would  be  reduced  to  M.  However,  another 
approach  can  be  taken,  giving  more  insight  into  the  problem  and  allowing 
the  use  of  a  general  procedure  for  any  delay.  Starting  from  the  result 
obtained  for  the  noncausal  filter,  observe  that  (C^  +  is  a  positive 
definite  matrix  and,  consequently,  it  can  be  uniquely  decomposed  at  [ll]  . 

C_  +  C  =  L  ♦  U 

— S  — —  — 

where  ^  is  a  lower  triangular  matrix  and  U  its  conjugate  transpose.  Thus, 

the  optimum  noncausal  filter  [G  .  1  assumes  the  form 
^  —A  one 

[g  1  =  c 

—A  ■'one  — S—  — 

where  U  ^  is  upper  triangular  and  L  ^  is  its  conjugate  transpose  or  simply 
transpose  if  F  is  a  real  process. 

Now,  consider  the  process 

W  = 

The  covariance  matrix  of^  is  given  by 

But 
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so  that 


C  =C+C  =LU 
— F  -S  — N  —  - 


=  i 


Thus  it  is  concluded  that  _L  performs  a  whitening  filter  operation  on 
Since  ^  is  a  lower  triangular,  this  is  equivalent  to  the  statement  that 
^  ^  represents  the  Gram-Schmidt  orthogonalization  process  on  the  data 
F.  The  problem  now  has  been  transferred  to  the  task  of  designing  an 
optimurxi  causal  discrete  filter  operating  on  the  white  noise  process  _W. 
This  can  be  done  by  first  examining  the  optimum  noncausal  filter  operating 
on  W,  which  is  given  by 


\gJ  L  =  C^U‘ 
—A-' one—  "  — 


The  product  C  U  can  be  decomposed  as  a  siim  of  a  lower  triangular  and 
an  upper  triangular  matrices 


CgU-'  =  [CsU-*\  +  [£sJJ'‘^u 


The  first  term  operates  only  on  past  and  present  values  of  W,  which  are 
uncorrelated  with  future  values  of  W,  This  means  that  the  best  that  can 
be  performed  when  no  future  values  of  W_  are  to  be  operated  upon,  is  to  take 
the  optimum  causal  filter  [  G^"l  as 
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Observe  that  the  procedure  of  taking  [G.  1  =  [C  U  L  1  is,  in 

general,  wrong,  since  ^  is  not  necessarily  white  (an  exception  occurs  when 

A  is  the  Karhunen-Loeve  transformation)  so  that  future  values  of  ^  are 

correlated  with  past  or  present  values  of  this  process. 

The  Suboptimal  Filter  with  Non- Zero  Delay.  When  a  delay  of  D  units 

is  allowed  in  the  operation,  it  is  possible  to  use  values  of  the  process  W, 

D  units  of  time  after  the  present  time.  Then  the  optimum  filter  under  the 

constraint  of  delay  D,  (G  .  )  has  the  expression 

A  oCl 


h  -  [C,U-  1 


— s—  -'LD 


(GA'od  = 


-s—  "LD 


where  [C  U  1  .  includes  the  elements  of  a  lower  triangular  matrix  and 
subdiagonals  up  to  the  delay  D. 

The  number  of  multiply  and  add  operations  involved  in  the  filter  matrix 
multiplication  ^  ^  when  a  delay  D(0  ^M-1)  is  allowed  can  now 

be  calculated.  The  first  (M-D-1)  rows  of  (from  top  to  bottom)  have 

(M- 1 ),,..,  (D+1)  nonzero  entries.  Each  of  the  last  (D+1)  rows  have  M 
nonzero  entries,  so  that  the  total  number  is  given  by 


(D+1)M  +  (D+M)  =  |(M(M+D)  +  (D+1)(M-D)) 


Calculation  of  the  Additional  Mean  Square  Error.  When  the  optimum 
noncausal  generalized  Wiener  filter  is  used,  it  can  be  shown  that  the  mean 


square  error  i 


is  [l]  . 
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=  t^A^onc  - 

The  fact  that  ^  )  is  a  positive  definite  matrix  guarantees  that  the 

o  In 

error  is  strictly  increased  for  any  nonmaxim\im  allowed  delay. 

The  mean  square  error  for  different  transformation  varying  the 
preassigned  delay  from  (1  to  15  has  been  calculated,  for  a  one-dimensional 
(16  X  1)  input  data  vector.  The  input  process  coTi'’ists  of  a  Markov  process 
signal  in  the  presence  of  white  noise  as  defined  in  Section  4.  1.  It  can  be 
observed  from  the  curves  of  Figure  1  that  only  a  small  price,  in  terms 
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Causal  Wiener  Filters 


of  mean  square  error,  is  paid  by  allowing  a  few  units  of  delay  or  even  no 
delay,  for  any  of  the  considered  transforms.  This  is  due  to  the  concentration 
of  energy  on  and  around  the  main  diagonal  of  the  filter  matrix.  The  most 
severe  increase  in  mean  square  error  occurs  for  the  identity  operator, 
followed  by  the  Hadamard  transform.  As  should  be  expected,  the 
Karhunen-Loeve  filter,  being  diagonal,  is  not  affected  by  the  restriction 
of  causality. 

However,  the  final  judgement  for  the  selection  of  the  best  transformation 
requires  a  joint  consideration  of  the  curves  of  Figure  1  and  the  number  of 
operations  involved.  Another  possible  element  for  consideration  comes 
from  a  detailed  analysis  (beyond  the  scope  of  the  present  paper)  of  numerical 
errors  in  computing  the  inverse  of  (C  +  ^  )  in  the  general  case,  for  differ- 
ent  transformations. 

It  should  be  observed  that,  although  the  present  method  results  in 
causal  operation  only  on  the  matrix  in  the  case  of  the  identity  operator, 
the  overall  transformation  is  causal. 

Finally,  consider  the  fact  that  the  _L  ^  factorization  method  (also  called 
square  root  factorization  or  Cholesky  method)  suggests  an  efficient  way  of 
inverting  the  covariance  matrix  (^  +  C  ),  namely,  obtaining  the  lower 

O  iN 

-1  .T  -1 

triangular  matrix  Ij,  inverting  it  and  multiplying  (L  )"'  by  Ij 

Experimental  Results.  The  suboptimal  filter  has  been  used  to  enhance 
a  256  by  256  pixel,  64  grey  level  picture  of  a  toy  tank  shown  in  Figure  2a. 

In  Figure  2b,  white  Gaussian  noise  has  been  added  under  the  condition  of 
unity  signal  to  noise  ratio.  This  was  done  by  setting  the  variance  of  the 


(c)  Optimal  Hadamard  Filter 
(delay  =  15,  SNR  =  1.0) 


Causal  Hadamard  Filter 
(delay  =  0,  SNR  =  1.0) 


(a)  Original  Picture 

(b)  Original  plus  Gaussian  Noise 
(SNR  =  1.0) 

Reproduced  from 
bes!  available  copy. 

noise  equal  to  the  measured  variance  of  the  pixel  values.  Figures  2c  and 
2d  show  the  comparison  of  the  performance  of  the  optimal  (delay  =  15)  and 
suboptimal  (delay  =  0)  filters,  with  Hadamard  preprocessing.  The  correlation 
coefficient  giving  the  best  result  was  found  to  be  equal  to  P  =  0.  5  for  this 
particvilar  picture.  Hardly  any  difference  between  the  two  pictures  was 
noticed. 

Conclusion.  The  concept  of  causal  filtering  has  been  extended  to 
generalized  Wiener  filter  computations.  The  method  is  based  on  a  factori¬ 
zation  of  the  covariance  matrix  of  the  signal  plus  noise.  It  has  been  shown 
that  only  a  small  price  in  mean  square  error  is  paid  even  when  small  or  iio 
delay  is  allowed  on  the  filter  matrix  multiplication.  The  comparison  of  the 
performance  of  the  optimal  and  suboptimal  filters  indicate  that  even  though 
the  number,  of  operators  for  the  filter  mulUplication  is  reduced  to  almost 
half  for  the  zero  delay  filter,  very  little,  if  any,  visual  difference  is  detected 
on  the  filtered  pictures  for  different  unitary  transforms. 
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4.  3  Homomorphic  Scalar  Wiener  Filtering  ^or  Image  Restoration 
Robert  H.  Wallis 

Images  degraded  by  the  effects  of  non-uniform  illumination  can  be 
dramatically  restored  by  using  a  non-linear  filtering  technique  which 
utilizes  the  response  of  the  human  visual  system  [l*]  .  Experimental 
evidence  indicates  that  the  human  eye  behaves  in  accordance  with  the  model 
shown  in  Figure  1.  This  model  explains  the  eye's  ability  to  accommodate 
an  enormous  dynamic  range  of  intensities.  Consider  the  input  to  the 
system  as 

I  =  i  •  r  (1) 

xy  xy  xy 

where  i  represents  the  illumination  of  the  scene  and  r  its  reflectivity, 
xy  xy 

The  logarithm  transforms  this  product  into  the  sum 

logd  )  =  log(i  )  +  log(r  )  (2) 

xy  xy  xy 

Since  the  illumination  of  natural  scenes  varies  gradually  across  the 

image  while  the  reflectivity  contains  much  sharp  detail  (edges,  etc.  ),  the 

high  pass  filter  tends  to  extract  the  scene  from  the  shadows.  By  treating 

the  log  (i  )  term  as  additive  inteference,  a  processor  which  imitates  the 
xy 

human  eye  can  be  constructed,  reducing  the  effects  of  nonuniform  lighting 
of  the  scene.  A  block  diagram  of  this  processor  is  provided  in  Figure  2. 

The  homomorphic  filter  adds  an  exponential  function  to  the  visual  model 
in  order  to  neutralize  the  logarithm.  The  log  function  transforms  inten¬ 
sities  into  densities.  The  filtered  densities  are  then  transformed  back 
into  intensities  for  viewing.  Since 
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Figure  4.  3-?.  Homomorphic  Filter 
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(3) 


logd'  )  = 
xy 


log(i'  )  +  log{r'  ) 
xy  xy 


can  be  treated  as  signal  plus  additive  interference,  the  linear  filter  which 
best  estimates  the  signal  in  a  mean  square  sense  can  be  found  using 
Wiener  filtering  theory. 

Wiener  filtering.  Wiener  filtering  is  a  technique  which  estimates  a 
signal  in  the  presence  of  additive  noise,  or  interference,  with  a  minimum 
mean  square  error.  Consider 


+  n  (4) 

where 

£  =  available  data  vector 
£  =  signal  vector  to  be  estimated 
n  =  interference  vector 


All  vectors  are  of  length  N.  For  simplicity  let  the  structure  of  the  estimator 
be  as  follows 


(5) 


That  is,  the  estimate  of  the  k  element  of  the  vector  £  is  given  by  the 
vector  dot  product 


h,y, 


(6) 


The  problem  then,  is  to  solve  for  the  h  that  minimizes  the  squared  error 


E(e^)  =  E[h"^£  -  £j^l^ 


(7) 
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Where  E  denotes  expected  value.  Defining  V,  to  be  the  gradient  with 

h 

respect  to  h,  it  is  required  that 


therefore 


V^{E[h\  -  Sj^ri  =  0 


T  T  T 

h  E{y_y_  )  =  E{s^  ) 


h  =  [E(^^^)]  ^E(s^) 


Assvirning  the  signal  to  be  uncorrelated  with  the  interference,  and  that 
both  are  zero  mean,  then 


(11a) 


E(y  y  )  ---  C  +  C 
s  n 


(lib) 


where  C  and  C  are  the  covariance  matrices  of  the  signal  and  interference 
8  n  “ 

respectively.  Therefore 


h  =  [C  +  C  r  E(s,  s) 

—  '•  e  Ir  — 


An  Example  of  Homomorphic  Wiener  Filtering.  For  simplicity  assume 
that  both  the  illumination  densities  (interference)  and  the  reflectivity  densities 
can  be  modeled  by  Markov  processes,  i.  e. , 


C  {i,j)  =  E{n.n.)  = 
n  IT  n  n 


(13a) 
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C(i,j)  =  E(s.=.)  = 


lk-!l 


lk-2l 


E(s 

“k 


5?^ 

-*  f 


2 

o 

s 


1 

P 

2 

P 


(13b) 


(13c) 


Since  the  interference  is  highly  correlated  (slowly  varying)  and  the  signal 
less  correlated,  p  should  be  chosen  larger  than  P  .  The  following  para¬ 
meters  were  selected  somewhat  arbitrarily. 


a 

=  0  =1 

(14a) 

n 

s 

P 

=  .9 

(14b) 

n 

P 

=  .6 

(14c) 

s 

k 

=  8 

(14d) 

N 

=  16 

(14e) 

Thus,  the  resulting  h  vector  is  16  elements  long  and  estimates  the  8th 
element  of  the  ^^-vector.  This  scalar  filter  can  be  used  to  filter  large 
records  of  data  by  convolving  h  with  a  long  data  vector.  This  can  be 
thought  of  as  multiplication  of  the  data  array  by  a  large  circulant  matrix. 


-72- 


filtered 

signal 


filter 

matrix 


data 


Each  row  of  the  filter  matrix  is  a  shifted  version  of  h  centered  at  the 
diagonal.  Of  course,  the  first  and  last  N/2  elements  of  the  £  wector  will 
be  inaccurate  estimates,  however,  for  the  large  arrays  used  in  image 
processing,  this  error  is  negligible. 

The  actual  convolution  is  carried  out  most  efficiently  using  the  Fast 
Fourier  Transform  (FFT)  algorithm  [s"].  The  processing  is  extended  to 
two  dimensions  by  using  a  two  dimensional  FFT  whose  frequency  domain 
filter  function  is  the  product  of  two  one-dimensional  filter  functions,  i.  e. , 

=  H(£^).H(£^)  (16) 

where  f^  and  f  are  spatial  frequencies  in  the  x  and  y  directions. 

Although  the  filter  parameters  have  been  chosen  somewhat  arbitrarily, 
the  results  are  quite  good.  Figures  3a  and  3b  present  original  and  restored 
images  of  a  tank  sampled  on  a  256  by  256  grid.  Note  that  the  restored 
version  shows  rruch  detail  not  visible  in  the  original,  and  seems  to  be 
more  uniformly  illuminated.  It  is  expected  that  more  careful  modelling 
will  yield  even  better  results. 
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(a)  Original 


(b)  Reconstruction 


Figure  4.  3-3.  Homomorphic  Filter 
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Conclusion.  This  is  a  preliminary  study  and  neglected  many  important 
considerations  such  as  the  calibration  of  the  intensities  in  the  original 
image.  Questions  such  as  "what  effect  does  the  length  of  the  filter's 
impulse  response  have  on  its  performance?"  remain  unanswered.  The 
results  however  are  quite  encouraging  and  future  study  on  extensions  of  the 
technique  are  planned. 
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4.4  Image  Restoration  By  Space-Variant  Decomposition 
Alexander  A.  Sawchuk 

In  the  well-known  linear  superposition  model  for  incoherent  imaging, 
the  degradation  from  object  to  image  is  given  by 

+® 

J'(x)  =  J  h(x,uK>(u)du  (1) 

^00 

where  (3(u)  is  the  original  object  intensity  function,  »P(x)  is  the  image 
intensity  recorded  by  the  system,  and  h(x,  u)  is  the  response  at  image 
coordinates  x  =  (x  ,  x  )  to  a  unit  impulse  at  u  =  (u  ,  u  )  in  the  object  coordin- 

1  ^  Xu 

ates.  In  general,  the  form  of  h(x,  u)  varies  with  the  position  of  the  object 
impulse  u  and  is  called  a  space-variant  point- spread  function  (SVPSF)  in 
the  context  of  image  processing.  If  h(x,u)  can  be  written  as  a  function  only 
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of  the  difference  x-u,  then  h(x,  u)  is  a  space-invariant  point-spread 
function  (SIPSF)  and  (1)  reduces  to  a  convolution. 

A  large  number  of  enhancement  and  restoration  methods  are  available 
to  process  degraded  images,  but  so  far  most  of  them  are  of  practical  use 
only  when  the  system  degradation  (1)  is  space -invariant.  Restoration  in 
the  general  space-variant  case  is  theoretically  possible  using  brute  force 
sampling  and  digitization  of  the  functions  in  (1)  followed  by  solution  of  the 
resulting  tensor  equation,  but  the  enormous  computational  capacity  required 
for  practical  realization  makes  such  a  technique  useless  in  most  cases. 

With  these  facts  in  mind,  it  is  evident  that  the  restoration  of  space-variant 
degradations  might  be  easily  accomplished  if  the  SVPSF  could  be  uniquely 
related  to  an  equivalent  space-invariant  system.  This  approach  to  image 
restoration  is  referred  to  as  space-variant  decomposition. 

Given  a  general  linear  SVPSF  describing  an  optical  system  degradation, 
the  basic  idea  behind  decomposition  is  to  express  the  SVPSF  as  an  equivalent 
cascade  of  invertible  geometrical  distortion  operations  with  space -invariant 
systems.  Figure  1  illustrates  the  concept  of  a  space-variant  system  h(x,u) 
and  its  equivalent.  Because  the  distortions  are  one-to-one  mappings  of 
points  between  coordinate  systems,  they  are  invertible  even  though  they 
may  be  nonlinear  in  the  spatial  coordinates.  The  center  block  is  a  space- 
invariant  system  hj(^-v)  which  gives  the  transformed  intermediate  image 
J  (£)  by  operating  on  the  transformed  object  O  (v).  Once  h(x,  u)  is  reduced 

Z  V 

to  this  form,  the  restoration  proceeds  ''-y  inverting  the  distortions  and 
applying  computationally  simple  space- invariant  estimation  or  inverse 
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filtering  to  recover  the  object  fiinction  [sl .  This  idea  is  summarized  in 
Figure  2  and  is  a  generalization  of  techniques  reported  previously  [2^,  [3^  . 

The  procedure  for  testing  whether  or  not  the  decomposition  is  possible 
for  a  given  SVPSF  involves  manipulating  h(x,  u)  into  a  functional  form  in 
which  the  geometrical  distortions  and  space -invariant  hj(2-v)  are  more 
obvious.  Simply  stated,  if  h(x,u)  is  expressible  in  the  form 

h(x,u)  =  a(x)P(u)hj(c(3:l  -  b(u))  (2) 


where  a(x),  3(u),  c(x)  and  b(u)  are  one-to-one,  continuous,  and  invertible 
functions,  with  a  (x)  ^  0  and  P  (u)  ^  0  over  their  domains,  and  hj(r)  is  a 
Fourier  transformable  function,  then  space-variant  decotnposition  and 
restoration  are  possible.  The  SVPSF  h(x,u)  may  also  be  put  in  alternative 
related  forms  from  which  restoration  follows  [  l].  Writing  h(x,u)  in  the 
form  (2)  immediately  gives  the  geometrical  transformations 


z  =  c(x) 


V  =  b(u) 


and  the  SIPSF  hj(r).  Here  the  intermediate  functions  are  given  by 


0  (v)  =  P(b  \v))J  i(x)0(b  \v)) 


JJ(c  (z)) 

(,)  .  - 

a(c  (z)) 


where  J,  _l(v)  is  a  Jacobian  function. 


Figure  4.4-2.  General  Space-Variant  Decomposition 


By  examining  physical  sources  of  optical  image  degradation  from  a 
system  point  of  view,  it  has  been  fotmd  that  the  decomposition  can  be 
derived  for  a  large  class  of  space-variant  imaging  systems  [l],  [2^,  [5*|  . 
The  application  of  this  technique  to  incoherent  motion  blur  is  discussed  in 
the  following  section,  while  the  remainder  of  this  section  is  concerned  with 
the  restoration  of  optical  aberrations. 

In  optical  imaging,  degradations  due  to  the  system  aberrations  are 
usually  of  greater  magnitude  than  those  due  to  diffraction.  An  inherent 
difficulty  in  attempting  restoration  arises  because  these  aberrations  are 
usually  quite  space-variant.  Robbins  and  Huang  [3*1  have  derived  the  form 
of  the  aberration  point- spread  function  as 


h(r.,e.;r  ,9  ) 

1  1  o  o 


-t 

r , 
o 


r.  cos(0,-0  )  -  r 

JL _ 1  O 

n 


r.  sin(0.-0  )' 
-+  1,  ^  ° 


n 


(6) 


where  (J'^f  ®  (r.,0^)  are  object  and  image  respectively,  expressed  in 

polar  coordinates  for  simplicity  due  to  the  circular  symmetry  inherent  in 
many  optical  systems.  In  this  formulation  the  function  h'‘'(x.,y^)  is  the 
response  to  an  input  impulse  at  standard  position  r^^  =1,  9^  =  0,  and  the 
expression  for  the  overall  system  operation  becomes 

2n  « 


Ar  .Q  ]  =  J  J  h(r  ,9, 
0  0 


r  ,  9  )C^  (r  ,  9  )r  dr  d9 
o  o  o  o  o  o  o 


(7) 


with  n  =  0  representing  spherical  aberration,  n  =  1  coma,  and  n  =  2  astigmatism 
and  .;ui  .  ature  of  field. 
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With  knowledge  of  the  degrading  expression  (6),  it  is  often  possible  to 
use  the  decomposition  technique  for  space-variant  restoration.  Robbins 
and  Huang  [s]  reported  a  method  for  the  restoraticn  of  coma  (n  =  1)  which 
is  essentially  a  special  case  of  the  general  technique.  By  making  the  changes 
of  variable  (the  geometrical  distortions  in  Figi  re  1) 


In(r.) 

1 

ln(r  ) 
o 


one  obtains 


h^  e 


zi-vi 


cos(z2-V2);  e 


..^l-^l 


sin{z^-v^) 


i)  =  V"i-''r 


where  the  1/r  factor  from  the  right  hand  side  of  (6)  is  absorbed  into 
o 

0^(r  ,  0  ).  It  can  be  shown  [ll  that  a  Mellin  transform  method  which 
0  0 

accomplishes  the  same  coma  restoration  is  really  the  same  decomposition 
technique. 

The  application  of  decomposition  restoration  to  spherical  aberration 
(n  =  0)  and  other  aberrations  (n  =  2)  is  not  quite  so  straightforward  but  is 
possible  for  certain  special  cases  of  h*(x.,y.)  and  is  a  currently  active 
research  area. 

The  decomposition  method  is  also  useful  when  the  degradation  is  a 
space-variant  blur  caused  by  tilt  of  the  image  plane.  For  this  case,  the 
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SVPSF  can  be  represented  approximately  by  [s'] 


^  '-T-j  • 

u,  '1  1  ' 


(10) 


When  h'''(u  ,  u  )  factors  into  functions  h'‘'(u,)h;r(u  )  it  is  called  separable, 
and  for  this  large  class  of  functions  space-variant  restoration  of  tilt  can 
be  accomplished  by  decomposition. 

The  restoration  methods  described  here  are  very  general  and  useful 
in  both  noisy  and  noise-free  systems.  Future  work  will  be  directed  at 
determining  other  areas  of  application  and  implementing  digital  simulations 
of  the  technique.  Several  papers  reporting  on  the  details  of  this  work  are  in 
preparation  and  portions  will  be  presented  at  the  Fall  Convention  of  the 
Optical  Society  of  America  [s'] . 
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4.  5  Restoration  of  Motion  Degraded  Images 
Alexander  A.  Sawchuk 

The  coordinate  transformation  restoration  technique  described  in  the 
previous  section  is  of  particular  value  in  the  restoration  of  ima.ges  degraded 
by  motion.  In  a  recent  paper  [l^  an  equivalent  linear  space-variant  system 
containing  all  the  motion  effects  has  been  derived  given  a  mechanical  des¬ 
cription  of  the  motion.  The  mechanical  description  has  the  general  parametric 
form 


Uf  =  gj(x^,X2;t)  =  gj(x;t) 

(la) 

^2  = 

(lb) 

which  uniquely  relates  any  object  point  u  =  (u^.u^)  to  the  location  x  =  (x^^x^) 
of  its  image  in  the  fixed  frame  as  a  function  of  time  for  the  exposure  interval 
[0,T]. 

Although  the  equivalent  space-variant  point- spread  function  (SVPSF) 
which  can  be  derived  from  (1)  is  rather  complicated  [l],  [2]  it  reduces  to 
some  simply  expressions  for  common  types  of  motion.  For  parallel  object 
and  image  intensity  function  planes  which  translate  during  exposure,  equation 
(1)  takes  the  form 


gj(x;t) 

=  X  -m,(t) 

1  1 

(2a) 

g2(x;t) 

(2b) 

where  only  the  coordinate  difference  x-u  is  a  time  function,  and  the  m.(t) 

-  1 

describes  the  planar  translation.  This  particular  form  of  motion  function 
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leads  to  an  equivalent  motion  point- spread  function  which  is  space-invariant, 
yielding  a  simple  solution  to  the  restoration  problem. 

More  complex  space-variant  motion  blur  occurs  when  object  and  image 
planes  are  not  parallel.  For  a  large  class  of  aerial  imaging  systems,  the 
motion  functions  (1)  can  be  expressed  as 


bi(U,,Uj,) 

(3a) 

*’z'“r“2’ 

=  -  m^d) 

(3b) 

where  the  mj^(t)  express  a  translation  between  transformed  coordinate 
planes.  Making  a  change  of  variable 


'^l 

(4a) 

=  b^CUj.U^) 

(4b) 

^1 

=  c,(x,,X2) 

(4c) 

'^2 

(4d) 

immediately  puts  (3)  into  space- invariant  form  and  simultaneously  identifies 
the  coordinate  distortions  for  space -variant  decomposition  described  in 
Figure  1  of  Section  4.4.  The  general  restoration  method  of  Figure  2  of 
Section  4.4  is  then  applied  directly  with  a  space -invariant  estimator  or 
inverse  filter.  The  space -invariant  equivalent  system  h  (w  ,w_)  is  given  by 

±  ±  Cl 

1 


?  r 
^-1  2 


[(mj(t))  f  (m^ft))  ] 


t=mj\wj) 

t=m-^W2) 


mj(0)  s  Wj  ^mj(T) 

m2(0)  ^  w^  ^ni2(T) 


elsewhere 
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(5) 


Coordinate  transformation  restoration  has  been  applied  in  a  computer 

simulation  of  space-variant  aerial  imaging.  Figure  la  is  the  ideal  lunar 

object  scene  taken  by  a  Ranger  spacecraft.  Figures  lb  and  Ic  respectively 

show  the  motion  degradation  when  the  object  is  imaged  by  side  and  forward 

oblique  tilted  cameras  with  constant  velocity  translation  at  fixed  altitude 

during  exposure.  Jn  Figure  lb,  the  camera  looks  sideways  to  the  left  at 

an  angle  of  20°  from  vertical.  The  movement  smears  the  image  as  a 

function  of  position,  verifying  the  space-variant  nature  of  the  system. 

Figure  Ic  shows  the  image  x-ecorded  by  a  forward  oblique  system  tilted 

forward  in  the  direction  of  movement  by  an  angle  20°  from  vertical.  As 

before,  the  blurring  length  and  direction  changes  with  position,  and  in 

contrast  to  the  side  oblique  case  there  is  blur  in  both  the  x^  and  x^  directions 

even  though  the  movement  is  along  the  axis.  For  these  two  space-variant 

systems,  it  turns  out  that  the  first  geometrical  transformation  v  =  b{u) 

Figure  1  of  Section  4.4  is  not  present,  so  restoration  proceeds  by  first 

using  the  appropriate  inverse  geometrical  distortion  z  =  c(x)  to  obtain  the 

intermediate  image  J>  {z)  which  appears  as  Figure  Id.  Using  an  inverse 

z 

FFT  filter  then  gives  the  final  restored  object  in  Figure  le. 

Another  application  of  coordinate  transformation  restoration  occurs 
when  there  is  relative  rotation  between  parallel  object  and  image  planes. 

The  equivalent  SVPSF  may  be  derived  in  this  case  [l],  where  it  is  shown 
that  the  polar  coordinate  transformations 
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(a)  Original  Lunar  Object 


(b)  Degradation  by  Moving 
Side  Oblique  System 


t 

^2 


(c) 


Degradation  by  Moving 
Forward  Oblique  System 


t 


2 


(d)  Intermediate  Image 


(e)  Final  Restored  Object 


Figure  4.5-1.  Space- Variant  Image  Motion  Restoration 


{6a) 


X,  =  r  cos  8 

1  z  z 

x_  =  r  sin  8  (6b) 

Z  z  z 

u,  =  r  cos  6  (6c) 

1  V  V 

u_  =  r  sin  8  (6d) 

Z  V  V 

used  for  the  geometrical  distortions  in  the  model  of  Figure  1  of  Section  4.4 
convert  the  rotation  into  an  eqxiivalent  space -invariant  hj(^-v)  for  restoration. 
Several  computer  simulations  to  verify  the  rotational  restoration  technique 
are  now  underway,  as  well  as  investigations  of  other  types  of  decomposable 
space-variant  motion. 
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4.6  Linear  and  Nonlinear  Interpolation  for  Image  Restoration 

Anil  K.  Jain 

In  many  parameter  estimation,  data  interpolation  and  other  similar 
situations,  one  has  a  knowledge  of  the  lower  and  upper  bovinds  of  the  esti¬ 
mates.  For  example,  in  the  processing  of  image  data,  obtained  by  scan¬ 
ning  a  black  and  white  picture,  it  is  known  the  data  should  always  be  positive 
(since  the  intensity  of  reflected  light  at  any  point  is  ^0).  The  estimates 
obtained  via  optimal  filtering  may  violate  these  bounds  due  to  inaccuracies 
in  modelling,  e.  g.  ,  when  a  two  dimensional  image  is  modelled  by  a  linear 
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Markov  process,  given  the  nc.isy  observations  and  the  bounds  on  the  esti¬ 
mates  themselves  (or  hard  constraints).  The  problem  is  formulated  by 
considering  the  maximum  likelihood  function,  and  the  presence  of  these 
bounds  (or  hard  constraints)  leads  to  a  nonlin'  ar  filter  [l]  . 

Consider  the  class  of  two-dimensional  images  represented  by  the  follow¬ 
ing  first  order  vector  Markov  proce^a  model  [z]. 


-\+l 

(1) 

II 

c  x^  +  ,  k  =  0,  1,...,N 

(2) 

where,  are  (N  X  1)  state  vectors,  are  (R  x  1)  observation  vectors, 
and  are  (N  X  1)  are  zero  mean  uncorrelated,  white  Gaussian  noise 
vectors  with, 

e[u^^]  =  K6(k-s) 

E[\]igl  =  L6(k-s) 

Now,  suppose  there  exists  a  priori  knowledge  of  the  lower  bound  on  the 
observations,  viz.  , 


V  ^  ^  k  =  0, 


N 


(3) 


where  ^is  an  (R  X  1)  vector.  The  problem  then  is  to  determine  the  optimal 
estimates  of  3^,  given  the  noisy  observations  and  the  constraints  so  that 
the  conditional  probability. 
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(4) 


T)  = 


max  fp(2Ev  5  iSi-  +  IL  ^  =  0»  •  •  • .  N)} 

^,k=0,N 


is  maximized.  The  solution  of  this  maximization  problem,  obtained  via 
the  calculus  of  variations  after  using  the  Bayes  rule  and  the  probability 
densities  of  and  u  ,  gives  the  following  nonlinear  equations  of  the  inter- 

iC  iC 

polator  [l]  (or  the  Maximum  Likelihood  Filter). 


in.,  =  Ax.  +BKB 


where 


-  p-^x 
-  ^ 

=  0 

i  * 

,  vise  1 

I-  2L  i 

1  -  z(xj 


and 


^-£x.  r  1 


(5) 

(6) 

(7) 

(8) 


(9) 


(10) 


It  can  be  easily  checked  that  for^  =  -®  (i.  e. ,  in  the  absence  of  any  con¬ 
straints),  equation  (6)  reduces  to  a  linear  equation  and  equations  (5)  and  (6) 
represents  the  standard  opilmal  linear  interpolator  equations.  These  linear 
set  of  equations  subject  to  conditions  (7)  and  (8)  are  easily  solved  by  using 
a  Riccati  Transformation  [l]  .  Figures  1  and  2  show  the  original  and  noisy 
image.  Figure  1  shows  the  optimal  linear  interpolation  on  the  two-dimensional 
32  X  32  noisy  image  with  S/N  ratio  of  7/3.  The  S/N  is  improved  by  roughly 
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20.  5db.  Let  u.  represent  the  linear  interpolation  estimate.  The  nonlinear 
interpolator  is  determined  by  first  linearizing  equation  (6)  about  the  linear 
estimate  u. .  The  linearized  equations  are  then  solved  via  a  Riccati 
Transformation  as  before.  Figure  2  shows  this  interpolation  estimate.  The 
major  improvements  in  the  picture  are  on  the  left  and  right  sides  of  the  white 
square  (about  1,  1  db).  The  improvement  in  the  picture  depends  on  how  often 
the  constraint  in  violated.  Better  approximation  to  the  nonlinear  filter  may 
be  obtained  by  taking  more  iterations.  For  a  picture  of  size  M  x  M,  the 
storage  requirement  for  interpolation  are  nM(M+l)  locations,  where  n  is  the 
order  of  the  vector  x  [2].  However,  an  estimator  or  one  step  predictor 
corresponding  to  the  nonlinear  filter  would  require  only  n  storage  locations, 
thereby,  making  on  line  filtering  feasible. 
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4.7  Pseudocolor  Image  Enhancement  Techniques 

Harry  C.  Andrews,  A.  G.  Tescher,  Richard  P.  Kruger 

The  article  commencing  on  the  following  page  appeared  in  the  July, 

1972  issue  of  the  IEEE  Spectrum  as  a  tutorial  paper  describing  digital 

image  processing.  It  is  included  in  this  report  because  of  the  relatively 

new  results  on  the  use  of  ;  seudocolor  which  appeared  on  the  latter  part  of 
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the  article.  The  article  lists  a  variety  of  techniques  on  the  use  of  color 
as  an  aid  to  viewing  monochrome  imagery,  and  illustrates  a  few  of  the 
methods  by  eolor  reproductions.  Applications  of  the  techniques  might 
range  from  biomedical  to  photo reconnai sane e  imagery,  the  success  being 
related  to  the  applicability  of  the  given  pseudocolor  technique  to  the 
viewing  objective  desired.  Unfortunately, this  is  more  of  a  user  oriented 
criterion,  and  one  in  which  theory  has  yet  been  quite  limited.  Future 
pursuits  of  the  technique  will  be  directed  along  the  lines  of  maximizing 
the  entropy  of  the  pseudocolor  displays  and  utilizing  just  noticeable 
difference  (jnd)  uniform  chromaticity  display  scales. 
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5.  Image  Detection  and  Measurement  Projects 

The  imagi  detection  and  measurement  projects  comprise  a 
mixture  of  tasks  concerned  with  techniques  and  tools  for  formatting 
images  and  extracting  information  from  pictures.  The  problems  con¬ 
sidered  under  this  classification  include  the  registration  of  pairs  of 
images,  the  measurement  of  image  features,  and  the  detection  of  ob¬ 
jects  within  pictures.  In  many  instances  the  solutions  to  the  problems 
are  very  much  unique  to  the  problem,  but  il.ere  are  also  universal  pro¬ 
blems.  The  latter  have  been  given  primary  consideration  in  the  study. 

The  first  report  describes  the  preliminary  effort  in  the  appli¬ 
cation  of  recursive  estimation  techniques  to  the  detection  of  objects 
within  pictures.  During  the  restoration  of  images  by  recursive  filter¬ 
ing,  one  of  the  by  products  of  the  estiiTiation  is  a  measure  of  the  likeli¬ 
hood  that  the  image  contains  an  object  of  interest  possessing  a  certain 
statistical  description.  This  information  can  be  utilized  with  a  minimal 
amount  of  effort  to  determine  the  presence  of  objects  in  the  presence 
of  substantial  amounts  of  image  noise. 

Another  research  topic  of  the  past  six  mcnthc  has  been  a  study 
of  edge  detection  techniques.  Several  methods  have  been  surveyed,  and 
one  of  them  -  called  a  bug  tracer  -  has  been  implemented  by  software. 
This  edge  detector  has  been  evaluated  on  a  sequence  of  poor  visual  quality 
radiographic  images.  The  performance  of  the  edge  detector  is  quite 
promising. 

5.  1  Detection  of  Objects  in  Noisy  Pictures 
Nassar  E.  Nahi 

Simple  recursive  procedures  have  been  developed  to  detect  the 
presence  of  objects  described  by  statistical  measures  in  pictures  con¬ 
taining  additive  noise.  Let  the  output  of  an  image  scanner  bt  denoted  by 
the  scalar  discrete  time  function  z{k)  =  y(k)  +  v(k)  observed  for  k=l,  .  .  .  , 

N  where  y(k)  and  v(k)  are  independent  scalar  stochastic  processes.  Sta¬ 
tistical  properties  of  v(k)  and  y(k)  are  assumed  known.  The  function  y(k) 
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Image  processing 
by  digital  computer 

Image  restoration  and  enhancement  can  be  implemented  by 
means  of  large-scale  digital  computers,  with  striking  results 

Harry  C.  Andrews  University  of  Southern  California 
A.  G.  Tescher  The  Aerospace  Corporation 
Richard  P.  Kruger  University  of  Southern  California 


Recently  d'sital  computers  have  been  brought  to  bear 
on  the  problem  of  processing  natural-quality  imagery, 
and  with  surprising  success.  Optical  processes,  communi¬ 
cation  systems,  and  automatic  image  detection  are  three 
areas  in  which  computer  techniques  are  being  applied. 
Perhaps  the  most  exciting  areas  of  research  are  those 
involving  image  restoration  and  enhancement.  The  most 
visually  dramatic  of  these  developments,  known  as  pseudo 
color,  is  a  technique  for  generating  color  from  black- 
and-white  images. 

Computers  and  pictorial  data 

High-speed  general-purpose  digital  computers  have 
reached  the  stage  of  development  where  it  is  becoming 
feasible  to  perform  mathematical  or  algorithmic  pro¬ 
cesses  on  images  of  natural  photographic  quality.  From 
an  analysis  viewpoint,  an  image  is  nothing  more  than  a 
two-dimensional  function  whose  variations  in  gray 
scale  describe  pictorial  information.  Because  man  has 
been  somewhat  successful  in  analyzing  one-dimensional 
signals  with  computational  techniques,  it  is  not  unreason¬ 
able  to  expect  that,  because  of  the  disproportionately 
large  role  that  the  eye  plays  in  human  experiences,  there 
will  undoubtedly  be  an  increasing  use  of  computers  for 
manipulation  of  such  natural  images.  Indeed,  within  the 
past  five  years  a  tremendous  amount  of  computer  pro¬ 
cessing  of  photographs  has  occurred,  with  facilities  spring¬ 
ing  up  around  the  globe  to  process  anything  from  aerial 
photographs  to  X  rays,  and  almost  everything  imaginable 
in  between.  However,  common  to  most  such  prt"esses 
is  the  fact  that  the  pictorial  information  is  most  often 
represented  as  a  brightness  function  of  two  space  vari¬ 
ables  (x,  y),  thus  allowing  two-dimensional  analysis 
techniques  to  be  utilized  In  fact,  the  pictorial  data  may 
not  be  imagery  at  all  but  range-Doppler  planes,  sonar- 
gram  displays,  or  voice  prints  In  any  event,  utilizing  the 
basis  of  two-dimensional  analysis,  it  is  possible  to  de¬ 


velop  universal  computational  techniques  for  a  large 
variety  of  image-processing  applications.  In  fact,  many 
two-dimensional  techniques  are  simply  trivial  extrapola¬ 
tions  of  one-dimensional  processes,  for  which  consider¬ 
able  analysis  is  already  available.  This  article,  then,  at¬ 
tempts  to  outline  a  few  methods  utilizing  digital  com¬ 
puters  in  the  manipulation  of  two-dimensional  data  for 
image-quality  applications. 

Optical  processes 

Long  before  large-scale  digital  computers  were  avail¬ 
able,  optical  devices  were  being  mathematically  analyzed 
as  to  their  capability  of  representing  an  object  of 
interest  mapped  onto  a  photographic  plate.  This  map¬ 
ping  process  can  be  viewed  as  a  system  that  was  found 
to  be  linear  in  intensity  for  incoherent  illumination  and 
linear  in  amplitude  and  phase  for  coherent  illumination. * 
Thus,  if  the  object  can  be  described  as  a  two-dimensional 
function  o{x,  y)  then  the  photographic  image  becomes 

f(x,  y)  =JJ  h{^,  7),  X,  y)  o(i,  v)  dr)  (1) 

where  ({,  rj)  are  the  spatial  coordinates  of  the  object  of 
interest  and  /»({,  rj,  x,  y)  is  a  function  describing  the  eflFect 
of  the  optical  devices  on  distorting  the  object  in  generat¬ 
ing  the  image.  The  latter  function  has  often  been  referred 
to  as  the  point-spread  function  or  impulse  response  of 
the  system,  because  the  image  would  equal  the  point- 


The  reader’s  attention  is  directed  to  the  July  1972 
Proceedings  of  the  IEEE,  a  special  issue  devoted  to 
digitai  picture  processing. 
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COLOR  REPRESENTATION  of  the  left  ventricle  of  a  human  heart,  with 
a  computer-detected  systolic  (blue)  and  diastolic  (red)  ventricular  outline. 
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fiiiKtion)  were  ihe  objeel  ('t' imerest.  If  the  l•lOl^I-^prea^.l 
fiinetion  is  ineariani  lo  shilTs  iii  lire  delia  l'iini.ni)n.  then 


/)(i.  v.  \  )  /llA  1 


i) 


(2) 


which  suggests  that  the'  iiiiage  is  realh  a  etuic I'hitinii  I'f 
the  nbieet  with  the  iioint-spreael  fiiiictH'ii 


/(.V,  1  ) 


If  - 


I]) /i(.\  i.  1  ali/ci/'t  (Ta) 


/  t.v,  \  )  i)(  \ .  1  )  •  /;( \ ,  i  ) 


t.'b) 


Such  s\  stems  are  referreci  m  as  liiii'ar  isoplanai  I'l  shift- 
iiivanant.  It  is  kiieiwii  that,  lor  these  systems.  I  oiirier 
transforms  play  a  maior  role  in  analysis,  ami 

I  {II.  (  I  ()i/(.  I )  Ihii.  I )  t  Ji 


wheie  uppercase  l\:nelions  imply  Ihe  transform  domain; 
that  is. 


//la.  I  ) 


L.(  .h- 


\ )  e\p  1  — /(a.v  -  t  i  ll  (/v  (Iv 


(5) 


I  hns.  for  linear-shirt-in\;iriant  systems,  in  theory  it 
shcnild  he  possible  to  reeo\er  the  obieet  from  the  image 
h\ 


(Hii.  i  ) 


I  \ii.  i  ) 
IHk.  !■) 


(6) 


and  inverse  I  ourier  transforming  I'r,  eciimtilenily .  by 
dei.'inohmg  the  image  by  Ihe  peiinl-spread  runelion, 
I  hi  digit.i!  computer  has  hi  n  used  as  a  device  for 
malheniatieally  inverting  lundoiiigl  the  degradation 
elfeel  resulting  from  the  point-spread  function  of  the 


FIGURE  1.  Fourier  transform  of  an  image.  A— Original 
(Surveyor  spacecraft).  B— Logarithm  of  the  magnitude 
of  the  Fourier  transform.  C— Low-pass-filtered  original. 
D— High-pass-filtered  original. 


optical  imaging  system. Indeed,  the  power  and  gen¬ 
erality  of  the  computer  allow  far  greater  flexibility  for 
image  manipulation  than  was  First  envisioned  from  the 
optical  viewpoint. 

Digital  implementation 

In  order  to  explain  how  computers  can  be  used  ef¬ 
fectively  in  image  processing,  a  brief  digression  to  dis¬ 
crete  mathematics  may  be  helpful.  Because  we  will 
eventually  be  interested  in  other  than  optically  generated 
images,  we  will  use  the  notation  that  our  original  function 
will  be  /(.V,  ,v).  Then  in  the  computer  the  brightness  or 
amplitudes  will  be  stored  as  an  N  X  N  array  of  numbers, 
where  we  have  limited  ourselves  *o  square  pictures  for 
ease  of  analysis.  The  infinite  integration  Fourier  trans¬ 
form  will  be  replaced  by  a  two-dimensional  Fourier 
series  expansion 


1  .v-i  v-i  i -2iri  1 

X]  E  /(V,  ,v)  exp -!  (h.v  +  (y)> 

=  =  (  N  I 

(7) 

which  will  be  our  discrete  analog  of  the  continuous  Four¬ 
ier  transform.  Figure  1  presents  an  original  image  (256  X 
256  with  64  gray  shades);  its  Fourier  transform,  given 
in  Eq.  (7);  and  a  low-pass  and  high-pass  version  of  the 
original.  The  filtered  images  have  been  obtained  by 
multiplying  F(m,  r)  with  a  filter  function  H(ii,  i)  such  that 


1  .V  -  1  .V  - 1 

/(-V.  .V)  =  E  E  f'O'-  '•)  W(M,  (•) 

^  u  =  0  r  =  () 


X  exp 


Ivi 

N 


(iix  -T  r.v) 


and,  for  Fig.  1,  H(u.  v)  is  a  zero-one  filter  retaining  or 
suppressing  the  spatial  frequencies  respectively. 

Whereas  the  Fourier  transform  finds  considerable  use 
in  modeling  and  inverting  optical  systems,  o't'er  mathe¬ 
matical  transforms  might  be  of  interest  for  coding  or 
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|,  FIGURE  2.  Two  digital  communiration  systems.  A— Dif> 
^  ferential  puisa^ode  modulation  system.  B— Transform 
I  image-coding  system. 


pattern-recognition  activities.  One  such  transform  that 
has  been  extensively  investigated  for  these  activities*-*  is 
a  subset  of  the  Hadamard  transforms  known  as  the  Walsh 
transform.’  *  The  advantage  of  the  Walsh  transform  is 
that  it  is  comprised  of  ±1  entries  and  thus  is  implement- 
able  by  simple  additions  and  subtractions.  In  notation 
equivalent  to  that  of  Eq.  (7),  we  see  that 


I  N-l  N-l 

F(u,  o)  =  -  E  S  fix,  .VX- 

^  z-0  y-0 


A— t 

2  UiZ{  +  til/i 

-iy-»  (9) 


where  xu  and  C(  are  the  binary  bits  representing  the 
spatial  and  transform  variables,  respectively,  and  N  = 

2«  9 

In  investigations  of  the  statistics  of  an  image,  second- 
order  moments  (variances)  become  particularly  rdevant. 
Assuming  separability  of  the  horizontal  and  vertical 
statistics  for  computational  simplicity,  let  [<^J  be  the 
covariance  matrix,  computed  or  modeled,  for  the  hori¬ 
zontal  or  x-dimension.  Then  it  is  of  interest  to  measure 
the  statistics  of  the  u-dimension  (the  transformed  x- 
variable)  in  the  various  transform  spaces.  In  general,  for 
the  iV  X  unitary  transform  [A], 


l<t>A]  =  [>l]to[>l] 


(10) 


For  the  Fourier  case. 


2mux 

[Ap]  =  exp  (11a) 

N 


for  the  Walsh  case, 
Mwi 


"S  u.xi  1 

=  (-ly-” 


and  for  the  Karhunan-Loeve  case, 
[Ak.l]  =  [EJ 

where 

=  [X]  = 


(llb) 

(llc) 

(12) 


in  transmitting  the  picture  across  a  digital  communica¬ 
tion  link.  As  an  example,  the  images  of  Fig.  1  require 
256  X  256  X  6  bits  for  representation  in  what  has  come 
to  be  known  as  a  PCM  (pulse  code  modulation)  mode 
(i.e.,  six  bits  for  each  picture  sample).  Yet  from  observing 
the  image  it  is  evident  that  there  is  considerable  redun¬ 
dancy  or  correlation  between  adjacent  samples  and  that 
it  might  make  more  sense  to  code  the  differences  in  ad¬ 
jacent  samples  rather  than  the  sample  values  themselves. 
This  technique,  or  various  modifications,  has  come  to 
be  known  as  DPCM  (differential  PCM)  and  is  the  fore¬ 
runner  of  some  of  the  methods  used  in  the  Bell  Labora¬ 
tories  Picturephone.  Initially,  it  might  seem  that  more 
bits  (seven  for  the  six-bit  PCM)  are  required  because 
we  must  be  able  to  code  both  positive  and  negative  dif¬ 
ferences.  However,  the  statistics  of  the  differences  are 
such  that  a  vast  majority  of  the  time  the  changes  are 
very  small,  so  coding  these  differences  with  short  code 
words  results  in  an  overall  bandwidth  reduction.  Good 
image-quality  reconstructions  have  been  achieved  with 
rates  equivalent  to  two  bits  per  picture  sample.  Figure 
2 A  is  the  block  diagram  of  such  a  system. 

In  transform  image  coding  (Fig.  2B)  the  mathematical 
transform  of  an  image  is  taken  in  two  dimensions — 
possibly  using  one  of  the  matrices  of  Eq.  (11)  twice, 
once  on  the  rows  (x-direction)  and  once  on  the  columns 
of  the  image  (y-direction).  The  transform  domain  is  then 
investigated  in  regard  to  its  coding  efficiency.  If  the 
Karhunen-Loeve  transform  is  used,  the  transform  sam¬ 
ples  are  uncorrelated  and  all  redundancy  has  been  re¬ 
moved.  Simply  coding  each  transform  sample  with  a 
code  word  proportional  to  the  expected  variance  of  that 
sample  will  result  in  a  quite  efficient  coding  scheme. 
Unfortunately,  computation  of  the  eigenvector  matrices 
[EJ  and  [Ey]  is  almost  impossible  for  large-array  images. 

Surprisingly  enough,  the  Fourier  and  Walsh  trans¬ 
forms  provide  adequate  domains  for  independent 
sample  coding  even  though  they  are  somewhat  inferior 
to  those  obtained  by  the  Karhunen-Loeve  technique.  '  ‘ 
In  addition,  they  have  a  rapid  implementation  algorithm 
that  allows  2iV*  log2  N  operations  instead  of  2A”  opera¬ 
tions,  which  becomes  a  considerable  savings  for  large- 
size  images. 


and  [X]  is  a  diagonal  m  .irix  comprised  of  the  eigenvalues 

0f[(^x]. 

Computer  image  coaittg 

The  objective  of  efficient  coding  of  an  image  is  to  re¬ 
move  the  redundancy  (usually  statistical)  from  the  image 


Computer  image  detection 

Because  of  the  successes  developed  in  the  field  of  de¬ 
tection  theory  for  one-dimensional  radar  returns,  it 
seems  natural  to  try  to  extend  the  analysis  to  two  dimen¬ 
sions  in  the  detection  of  objects  in  an  image  or  field  of 
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view.  The  computer  c:in  again  be  enlisted  as  an  aid  to  this 
task.  It  i-  known  that  the  matched  lilter  maximi.tes  signal- 
to-noise  ra'ios  for  detection  purposes  I  he  lilter  function 
utilized  in  F.i[  (S)  then  becomes'- 


//(».  r) 


F*{ii.  i  ) 


e\p  -/(//i  -  tr/) 


for  an  image  /(.v.  v)  loetited  at  point  (;.  i?)  in  noise  with 
a  Fourier  power  spectrum  r).  Indeed,  the  Fourier- 


transform  plane  becomes  the  natural  one  for  implement¬ 
ing  the  matched  lilter.  since  it  is  the  onl>  domain  that 
provides  an  invariance  to  position  of  the  image  one 
sought.  The  lilter  has  the  property  that  it  will  detect  the 
image  anywhere  in  the  plane  or  field  of  \iew  but  it  will 
reject  rotatio.ns  of  that  image.  Fhiis.  it  is  possible  to 
discriminate  between  the  p  and  </.  the  two  being  ISd- 
degree  rotations  of  one  another  Figure  presents  the 
results  of  such  an  experiment,  where  the  image  of  in- 
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FIGURE  3.  Image  detection.  A— Signal.  D— Array.  C— Array  and  low-frequency  noise.  D— 
Array  and  high-trequency  noise.  E— Filter  output  after  thresuold  detection.  F— Perspective 
of  correlation  peaks. 


FIGURE  4.  Motion-detection  images.  A— Cineangiogram  of  a  left  ventricle  at  diastole.  B— 
Cineangiogram  of  a  left  ventricle  at  systole.  C— Cineangiogram  of  a  left  ventricle  showing 
inserted  catheter  before  dye  injection. 


fields  of  view  in- 

ftrpT  *  ^  variety  of  noise  situations.  In  all 

three  cases  the  three  (fs  are  successfully  detected  with 

»f  %»1  P»ks  dLplaSt  Rg 
Another  example  of  image  detection  is  that  of  motion 

heSn  E  '"““i”  P“'POSes  is  describS 

herein.  For  several  years  cardiologists  have  routinely 

recorded.and  manually  measured  c  vdiac  dynamics  using 
cineangiography.  ”-15  Principal  di¬ 
agnostic  attention  has  been  concentrated  on  extracting 

r"'""* ""'““I" 

since  this  chamber  supplies  the  systemic  circulation  A 
radiopaque  blood-soluble  contrast  medium  is  injected 
into  the  heart  chamber  using  a  catheter  inserted  into  a 
peripheral  artery  or  vein.  As  the  contrast  medium  fills  the 

of  thTh™  the  dynamic  motion 

fiL  rnmV  I  "  fl“o™«:ope/image-intensi- 

her  combination,  the  output  of  which  is  interfaced  to  a 

vidicon  or  orthicon  television  camera  and  monitor 
Quite  often  the  output  of  the  image  intensifier  is  simul¬ 
taneously  recorded  on  a  16-  or  35-mm  motion  picture 
a  a  rate  exceeding  30  frames  per  second.  This  film 
record  is  then  used  to  make  manual  estimates  for  a  curve 
of  chamber  volume  vs.  time.”-”  Simultaneously,  a  re¬ 
cording  of  patient  electrocardiogram  and  chamber  pres- 

ZJh  “  ‘™“'  “■'“-^vnctodo^  J™. 

tion  of  diagnostic  importance. 

The  most  tedious  aspect  of  this  analysis  is  the  detection 

n  racing  of  the  chamber  boundary  as  outlined  by  the 

injected  contrast  medium.  Several  research  approaches 


taken  show  promise  of  an  eventual  solution.  However 
most  of^em  still  require  varying  degrees  of  human  inter¬ 
action  Some  recent  work  suggests  that  it  may  i 

nlTf  ‘fi**  boundary  detection  com- 

Plaely.”  thus  permitting  rapid  calculatim  of  heretofore 
difficult-to-acquire  diagnostic  information  such  as  ores 
sure-volume  loops  and  ventricular-mass  estimates 
Figure  4  shows  several  digitized  cinefilm  frames  of  a 
left  ventricle  in  the  right  anteri  -  oblique  position.  An 
example  of  completely  automated  ventricular  boundary 
tracing  using  color  displays  appears  on  page  21.  Volume 
estimates  obtained  from  the  outlined  end  systolic  and 
diastolic  ventricular  cavity  shown  are  less  desirable  than 
a  complete  volume-time  curve.  However  these  two 
volume  estimates  are  often  the  only  diagnostic  informa¬ 
tion  currently  extracted  manually  from  the  film.  From 
these  two  outlines  the  end  systolic  and  end  diastolic  vol¬ 
umes,  ejection  fraction,  and  stroke  volume  are  easily 
estimated.”-”  -  easily 

Computer  image  restoration 

Thus  far  in  this  article  we  haye  been  discussing  various 
techniques  for  processing  images  for  coding  and  detec¬ 
tion.  However,  possibly  the  most  exciting  aspects  of 
computer  image  processing  are  those  of  image  restora- 
tion.  Traditionally,  the  computer  has  been  used  to  “re¬ 
store  an  image  to  its  original  quality  according  to  some 
mathematical  manipulation  that  is  intended  to  invert 
some  physical  degradation  phenomenon  experienced  in 
the  formation  of  the  image  in  the  first  place.  Thus  the 
concept  of  image  restoration  is  one  of  attempting  to  re¬ 
cover  losses  suffered  in  the  optical  imaging  system  due 
to  a  variety  of  degradations,  which  can  be  modeled-from 
the  very  simple  to  the  extremely  complex 
Unear-^ft.jnv,  riant  phenomena.  Probably  the  sim- 
plest  matheiTOtical  model  of  the  imaging  process  was 
that  given  eariier  where  the  image /(x,  j.)  is  the  output  of 
a  linear  system  whose  input  was  an  object  o(x,  y).  The 
object  suffers  degradation  introduced  by  a  systems  re¬ 
sponse  A(x,  y,  f,  ,)  and  for  isoplanar  systems  the  image  is 
the  two-dimensional  convolution  of  the  object  and  a 
point-spread  function.  Ideally,  it  is  desirable  that  the 
point-spread  function  be  a  Dirac  delta  function,  in  which 
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S” with 


f(x,  y)  =  i(x,  y)  r(x,  y)  log  f(x.  y)  =  log  i(x,  y)  +  log  itx,  y)  H{u.  v)  F  [leg  f(x.  y)]  h(x.  y)  •  log  f(x,  y) 


F[log  f(x.  y)]  =  F[log  i(x.  y)]  +  F[log  r(x,  y)]  f(x,  y)  =  exp  [h{x.  y)  •  log  f(x,  y)] 

=  exp  [h(x.  y)  •  log  iC-  y)  +  h(x.  y)  •  log  r(x.  y)] 

FIGURE  6.  Homomorphic  filtering  techniquo. 
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case  the  image  equals  the  object.  In  fact,  such  is  not  the 
case  and  the  power  of  linear  systems  theory  can  be 
brought  to  bear  on  the  problem  of  restoring  the  image 
to  as  close  a  representation  of  the  object  as  possible. 
Thus,  for  isoplanar  degrading  phenomena  such  as  de¬ 
focus,  certain  linear  motion,  simple  spherical  aberrations, 
certain  atmospheric  degradation,  and  additive  noise  the 
mathematical  tools  are  available  for  computer  correction 
of  such  processes.  Inverse  and  linear  least-square  filters 
are  typically  employed  for  image  restoration;  gradients, 
Laplacians,  and  higher-order  spatial  derivatives  are 
often  useful  approximants.  Finally,  it  is  possible,  within 
the  linear  theory,  to  investigate  resolution  beyond  the 
diffraction  limit  by  using  the  tools  of  analytic  continua¬ 
tion  and  extrapolation.  Figure  5  displays  an  example 
of  inverse  filtering  according  to  Eq.  (6),  where  //(«,  c)  was 
modeled  as 

//(«,  p)  =  (14) 

indicating  that  the  point-spread  function  tended  to  roll 
off  in  spatial  frequency  according  to  a  two-dimensional 
Gaussian  surface.  The  two  originals  (Figs.  5A  and  5C) 
are  digitized  X  rays  of  normal  and  abnormal  patients, 
the  latter  exhibiting  a  <everi  case  of  pneumoconiosis 
(black  lung  disease).”  The  enhanced  images  make  this 
fact  a  bit  more  evident.  The  same  images  are  then  used 
in  a  pseudo-color  experiment  to  be  described  in  the  next 
section. 

Linear-shift-varying  phenomena.  For  more  complicated 
processes,  the  anisoplanar  model  may  be  more  relevant, 
in  which  case  the  point-spread  function  h{x,  y,  rj)  is 
different  at  each  point  in  the  two-dimensional  plane. 
Consequently,  the  convolution  equation  is  invalid  and 
we  must  revert  to  Eq.  (1)  for  our  model  and  our  intcgral- 
cquaticT  (kffiUion  PhLnotrera  dt'scribablc  utAt  thev- 
assumptions  include  geometric-distortion  correction, 
nonuniform  and  turbulent  atmosphere;  nonlinear-motion 
degradation,  coma,  tilt  in  cylinAical  lenses,  and  many 
others.  Often  it  is  possible  to  find  nonlinear  transforma¬ 
tions  that  map  the  image  into  a  space  in  which  linear 
systems  theory  can  then  be  utilized  to  tortett  fot  the 
degradation,  allowing  an  inverse  nonlinear  transforma¬ 
tion  to  describe  the  restored  image.*" 

Ncnlinear  model.  Possibly  a  more  relevant  model  for 
certain  imaging  situations  is  that  of  a  multiplicative 
process,  in  which  the  image  is  assumed  to  be  comprised 
of  the  product  of  an  illumination  function  and  a  reflected 
function.  The  power  of  homomorphic  filtering  can  then 
be  enlisted  to  separate  the  processes  end  filter  multi¬ 
plicative  noise  in  the  traditional  linear-systems  technique. 

If  the  image,  Ax,y),  can  be  considered  to  be  the 
product  of  an  illumination  term  and  a  reflected  term,  then 


f(x.,y)  =  i(x,y)rix,y)  (15) 

and  if  the  spatial  frequency  response  of  the  illumination 
term  is  lower  than  that  of  the  reflected  term,  then  linear 
filter  theory  can  be  utifi’"!  f'  •■emove  the  illumination 
degradation.  to  separate  out  the  multiplicative 

process,  a  logarithmic  operator  can  be  applied  followed 
by  the  linear  filter  for  linear-shift-invariant  processes, 
Eq.  (6),  followed  by  an  exponentiator;  see  Fig.  6.  Signifi¬ 
cant  restorations  have  been  accomplished  with  this 
technique.** 

Confiputer  image  enhancement 

The  computer  can  be  used  to  “improve"  the  quality  of 
an  image  without  recourse  to  knowledge  of  degrading 
phenomena  resulting  in  the  concept  of  image  enhance¬ 
ment.  Both  mathematical  and  heuristic  techniques  are 
utilized  in  this  process,  with  major  emphasis  upon  the 
“viewing”  of  the  image  for  extraction  of  information 
that  may  not  have  been  so  readily  apparent  in  the  pre¬ 
enhanced  original. 

Psychophysics  of  vision.  Although  the  computer  can  be 
used  quite  successfully  to  “restore”  images  due  to  some 
correctly  modeled  degrading  phenomenon,  its  power 
does  not  stop  here.  Specifically,  it  is  possible  to  model  the 
human  viewing  process  in  such  a  way  that  an  optimum 
display  is  computed  for  the  human  as  a  communication 
sink.**  Work  is  still  somewhat  limited  in  this  area,  but 
certain  combinations  of  logarithmic,  linear,  and  ex¬ 
ponential  processes  seem  to  provide  more  pleasing  view¬ 
ing  than  would  normally  be  expected.  Indeed,  if  a  \’iable 
fidelity  criterion  could  be  developed  for  the  human  eye, 
then  all  the  powers  of  information  and  communication 
theory  could  be  directed  toward  optimal  use  of  that 
chisniul  hk*ifcucti  lb*.  cciTfiuUt  trrjqa^.'  aad  nn*n  TTih 
topic  has  interesting  biological  as  well  as  engineering 
ramifications,  for  if  we  could  mathematically  understand 
the  human  viewing  mechanism  we  might  then  add  in¬ 
sight  into  the  biological  processes,  and  vice  versa.  Given 
a  proper  fidelity  criter-  m  for  human  viewing,  “optimum 
image  coding”  coiAd  be  envisioned  in  ■which  image  data 
not  relevant  to  the  viewing  process  could  be  eliminated 
before  transmission,  thereby  saving  channel  bandwidth. 

Heuristics.  Heuristic  techniques,  for  lack  of  a  fvetter 
description,  for  improving  the  quality  of  images  have 
been  developed  with  some  success.  Unsharp  masking,” 
contrast  enhancement,  and  enspening  have  all  been 
utilized  in  computer  image  processing.  Most  such  tech- 
riquev  xurt  limzr  and  shift -invcrinni  and.  ihacfure,  fit 
well  into  the  linear-filter  theory  but  are  not  necessarily 
designed  to  correct  for  a  specific  degradation.  Simply 
stated,  the  eye  appears  to  respond  well  to  high  spatial 
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FIGURE  7.  Demonstration  of  the  histo¬ 
gram-equalization  technique,  which  can 
be  employed  to  aid  in  the  viewing  of  com¬ 
puter-processed  images.  A— Original 
histogram  (a  probability  density  func¬ 
tion).  B— The  same  histogram  after 
application  of  the  equalization  tech¬ 
nique.  C— Original  X-ray  photo  of  a 
dog’s  heart.  D— Equalized  image  of  the 
dog’s-heart  photo.  E— Original  illustra¬ 
tion  of  a  tank,  photographed  through 
cheesecloth.  F— Equalized  image  of 
tank-and-cheesecloth  photo.  Note  that 
the  equalized  images  shown  here  pro¬ 
vide  more  content  than  the  originais. 


FIGURE  8.  Nonlinear-transform  filtering.  A— Block  diagram  of  the  process.  B— 
Originai  tank  image.  C— Defocused  tank  image.  D— Defocused  tank  image  after 
processing  (Fourier  transform  with  a  =  0.5).  E— Defocused  tank  image  after  pro¬ 
cessing  (Fourier  transform  with  a  —  0.7). 


F(u,  v)  =  log  |F(u. 

Ffu.  =  IFfu.  or  F(u.  v)  =  lF(u.  v)l“e 
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FIGURE  9.  Schematic  diagram  for  the  generatton  of 
pseudo-cclor  images. 


frequencies,  which  is  what  these  techniques  tend  to  em¬ 
phasize.  In  addition,  some  nonlinear  averaging  combined 
with  thresholding^*  can  be  utilized  in  an  anisoplanar 
filter  for  local  region  correction  and  variation.  One  par¬ 
ticularly  exciting  technique  that  aids  th**  human  in  view¬ 
ing  computer-processed  images  is  the  use  of  his  :ogram 
equalization;  see  Figs.  7A  and  B.  Here  it  is  argued  that 
as  many  shades  of  gray  as  possible  should  be  used  in  the 
display  of  an  image.  Such  a  criterion,  equivalent  to 
maximizing  the  zero-order  entropy  of  the  output  display, 
indicates  that  the  gray  shades  would  be  unifoimly  dis¬ 
tributed.  The  technique  is  equivalent  to  a  nonlinear  op¬ 
erator  on  the  intensity  scale,  which  is  position-invariant. 
Figon-S  fC  F  present  examples  of  ajuvpnter  processing 
of  images  according  to  the  scheme.  In  some  sense,  it  is 
evident  that  the  equalized  images  provide  more  content 
than  the  originals.  Equalization  is  particularly  im¬ 
pressive  with  X  rays  or  images  with  heavily  biased  histo¬ 
grams  toward  one  end  or  the  other  of  the  gray-scale 
range.**  The  technique  can  be  interpreted  as  making 
subtle  changes  more  evident  in  the  regions  for  which 
such  changes  occur  most  frequently,  while  losing  subtle 
intensity  changes  m  less  frequently  occurring  gray  scale 
regions. 

A  second  technique,  which  has  become  somewhat 
useful  for  investigating  sniaH  variations  in  usually  high- 
frequency  regions  of  an  image,  is  to  employ  what  could 
be  considered  the  inverse  of  the  homomorphic  filter 


FIGURE  10.  Chromaticity  triangles.  Left:  Maximum 
energy  constraint.  Right:  Maximum  component  con¬ 
straint. 


shown  in  Fig.  6.  Indeed,  it  has  been  found  that  a  nonlinear 
filter  in  the  transform  domain  that  tends  to  suppress 
large-valued  terms  (thereby  nonlinearly  enhancing  small¬ 
valued  terms)  results  in  rather  amazing  image  recon¬ 
structions.  Because  most  of  the  energy  of  an  image  is 
concentrated  in  the  lower-frequency  portion  of  the  Four¬ 
ier  domain  or  in  the  lower-sequency  portion  of  the 
Walsh  domain,**  if  a  logarithmic  nonlinear  device  is 
used  the  smaller  values  are  less  significantly  reduced 
than  the  larger.  These  small-valued  coefficients  occur  in 
the  higher-frequency  or  higher-sequency  regions.  To  be 
sure,  noise  and  weak  signals  are  both  enhanced  at  the 
expense  of  reducing  the  dominant  portion  of  the  lower 
and  more  uninteresting  frequencies.  The  examples  of 
Fig.  8  illustrate  the  point  of  nonlinear  suppression  in  the 
transform  domain.  Note,  carefully,  that  although  the 
image  quality  has  been  changed,  it  is  now  possible  to  view 
details  that  were  not  available  in  the  original  or  in  the 
defocused  version.  (The  wheels  of  the  tank  are  particu¬ 
larly  significant.) 

The  role  of  color  in  digital  image  processing  and  pseudo 
color.  Virtually  all  research  in  digital  image  processing 
to  date  has  involved  only  intensity  or  black-and-white 
images.  The  human  eye,  however,  perceives  color  as  well 
as  bri^ness  infotmation,  and  it  is  evident  that  the 
processing  techniques  should  eventually  include  color 
images  as  well.  We  can  envision  many  fields  where  color 
information  may  be  of  particular  value;  however,  we  will 
restrict  our  attention  to  a  specific  area— which,  we  be¬ 
lieve,  has  significant  benefits.  This  processing  technique 
involves  the  use  of  color  in  the  presentation  of  intensity 
images;  hence  the  name  “pseudo  color.”  The  pseudo¬ 
color  experiments  described  here  were  conducted  re¬ 
cently  at  the  University  of  Southern  California's  Im¬ 
age  Processir.g  Laboratory. 

It  is  well  knv.'vn  that  most  spectral  colors  can  be  repro¬ 
duced  to  the  human  viewer  by  an  ad»httve  color  system 
consisting  of  three  primary  colors  such  as  red,  green,  and 
blue.  Consequently,  for  the  appropriate  color  computer 
processing,  the  image  is  represented  by  three  irttensity 
matrices,  >•),  /©(x,  y),  and  Ib{x,  y),  where  x,  y  are 
the  spatial  locations  and  the  subscripts  R,  G,  and  B  de- 
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FIGURE  11.  “Simple”  pseudo  color  for  biomedical  images. 
Left:  Normal  lung  enhanced.  Right:  Lung  with  pneu¬ 
moconiosis,  enhanced. 


Reproduced  from 
best  available  copy. 


FIGURE  12.  "Conventional”  pseudo-color  schematic  dia¬ 
gram. 


note  red,  green,  and  biue  primary  components  respec¬ 
tively. 

Let  us  assume  that  each  primary  component  can  be 
represented  by  64  intensity  levels.  (The  hardware  device 
used  to  display  the  color  images  corresponds  to  this 
assumption.)  Clearly,  each  pixel  (picture  element)  can 
assume  64^  distinct  values,,  corresponding  to  all  color 
and  intensity  combinations.  Consequently,  expanding  this 
argument,  the  amount  of  information  that  an  image  can 
convey  to  the  observer  can  be  expanded  64'“-  or  4096- 
fold  using  the  indicated  hardware  device  as  compared 
with  similarly  quantized  black-and-white  images. 

Does  the  human  eye  accept  this  simple  mathematical 
argument?  The  fortunate  fact  that  the  answer  is  affirma¬ 
tive  provides  the  basic  motivation  for  research  in  pseudo 
colors.  Numerous  experiments  have  indicated  that  al¬ 
though  the  human  eye  can  differentiate  only  perhaps 


FIGURE  13.  “Conventional”  pseudo  color  for  biomedical 
Images.  Left:  Original.  Right;  Histogram  equalized. 
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FIGURE  14.  Processing  for  spatial-frequency  pseudo-coloi  ing, 


FIGURE  15.  "Unconventional”  pseudo-color  generation  for  tank  image.  Upper  left:  Com¬ 
posite  image.  Upper  right:  Green  n-imary.  Lower  left:  Blue  primary.  Lower  right:  Red 
primary. 
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FIGURE  16.  Spatial-frequency  content  of  the  tank-and-cheesecloth  image.  Left:  Composite 
image  with  red  content  removed.  Right:  Composite  image. 
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two  or  three  dozen  brightness  levels  in  an  image,  it  is 
able  to  separate  thousands  of  various  colors.’’ 

The  schematic  diagram  of  pseudo-color  generation  is 
shown  in  Fig.  9.  The  “process  and  decompose”  step 
essentially  includes  a  digital  computer,  where  the  ap¬ 
propriate  algorithms  are  implemented.  To  reiterate,  the 
objective  is  to  aid  human  interpretability  by  representing 
a  black-and-white  image  in  color  and  hence  take  ad¬ 
vantage  cf  the  inherem  capabf'ity  of  the  human  visual 
system  to  absorb  color  information.  This  technique 
actually  generates  color  where  none  exists. 

As  an  example  of  the  range  of  color  variation  that  is 
possible  by  our  display  device  (essentially  a  digitally  con¬ 
trolled  color  television  monitor),  the  two  chromaticity 
triangles  of  Fig.  10  were  generated.  At  the  left  all 
possible  colors  are  shown,  subject  to  the  energy  con¬ 
straint  that  Jhe  triplet  of  values  for  each  pixel,  repre¬ 
senting  primary  components,  add  up  to  a  constant. 
This  constraint  is  changed  so  that  without  altering  the 
ratio  of  primary  component  values  we  scale  each  pixel 
in  such  a  way  that  the  largest  of  the  primary  components 
has  the  same  value  for  each  pixel,  as  shown  at  the  right. 
Thus  the  second  triangle  contains  ail  possible  colors  at 
the  maximum  brightness  level  consistent  with  the  capa¬ 
bility  of  the  display  device.  Ti...  y...tices  of  the  two  tri¬ 
angles,  which  are  identical,  represent  the  three  primary 
colors.  Despite  the  severe  energy  constraint  imposed 
on  these  chromaticity  triangles,  more  than  2000  distinct 
pixel  values  are  contained  in  each  triangle.  In  contrast, 
the  black-and-white  images  discussed  previously  consist 
of  only  64  levels  of  gray. 

Historically,  the  earliest  implementation  of  pseudo 
color  can  be  found  in  cartography.  We  rarely  appreciate 
how  difficult,  if  not  impossible,  it  would  be  to  create  a 
complex  map  without  the  use  of  color.  Computer  im¬ 
plementation  of  pseudo  color  is  limited  only  by  the 
capabilities  of  the  display  device  and  the  imagination  of 
the  researcher.  The  approaches  may  vary  from  the  very 
simple  to  those  for  which  computer  implementation  is 


mandatory.  One  very  simple  technique,  which  does’  not 
mix  the  primary  colors,  involves  mapping  a  separate 
portion  of  the  gray  scale  into  maximum  brightness  for 
each  primary  component.  Figure  11  shows  examples 
for  this  case.  It  is  interesting  to  compare  these  X-ray 
images  of  a  normal  and  a  diseased  lung  with  their  black- 
and-white  counterparts  (Figs.  5B  and  5D). 

A  somewhat  more  complex  pseudo-color  gen-ration 
technique,  still  in  the  realm  of  the  conventional  approach, 
is  indicated  in  Fig.  12.  Mixing  of  the  primary  components 
in  this  case  is  performed  according  to  the  indicated  map¬ 
ping.  Here,  the  pure  primaries  occur  only  at  the  extreme 
ends  and  exact  center  of  the  gray  scale.  Figure  1 3  shows 
the  appropriate  examples  and  should  be  compared  with 
Figs.  7C  and  7D. 

If  maximum  entropy  is  the  objective,  as  in  the  case  of 
histogram  equalization  of  black-and-white  images,  a 
rather  complex  experiment  must  be  performed.’*  A  step 
in  this  direction  is  indicated  in  Fig.  14. 1'he  final  display 
of  this  experiment  (see  Fig.  8B)  is  shown  in  Fig.  15.  Note 
the  emphasis  of  the  wheels  and  the  separation  of  the  in¬ 
focus  foreground  from  the  blurred  background.  As  a 
final  example.  Fig.  16  utilizes  the  same  pseudo-color 
technique  to  demonstrate  the  spatial-frequency  content 
of  the  image  of  Fig.  7E. 

Clearly,  a  great  deal  more  research  is  required  before 
the  potentials  of  pseudo-color  displays  can  be  realized 
fully.  More  hardwt  re  development  is  also  required.  How¬ 
ever,  the  superior  i.'apr.bility  of  the  human  visual  system 
to  perceive  color  differences  provides  a  very  strong  in¬ 
centive  to  continue  research  in  this  field. 

in  conclusion 

It  should  be  emphasized  that  only  a  few  examples  of 
computer  techniques  in  image  processing  have  been  pre¬ 
sented  in  this  article.  Many  other  existing  applications 
have  been  excluded  either  for  space  reasons  or  because 
of  the  authors’  technical  limitations  in  certain  fields. 
However,  it  should  be  apparent  that,  with  powerful 


fL 


Andrews,  Tescher,  Kruger — Image  processing  by  digital  processing 


31 


digital  computers,  more  and  more  image-procesdng 
applications  will  arise  in  areas  yet  untouched  by  the 
techniques  covered  here.  ** 
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is  called  the  signal  process  and  v(k)  the  noise  process.  Given  z(k), 
k=l, .  .  . ,  N,  the  problem  is  to  decide  on  the  presence  or  absence  of 
the  signal. 

In  formulating  the  detection  problem  let  be  the  hypothesis 
that  the  signal  is  present,  and  be  the  hypothesis  that  the  signal  is 
absent. 


Hj  :  z(k)=y(k)  +  v(k) 


H  :  z(k)  =  v(k) 
o 


K=1,...,N 


Define  f.,  j=0,  1  as  the  joint  probability  density  of  z(l),  .  .  . ,  z(N)  when 

H.  is  true  and  define  the  corresponding  likelihood  ratio  as  t(N)=ln  — — 

^  o 

If  t(N)  exceeds  some  appropriate  threshold,  is  accepted;  other¬ 
wise  Hq  is  accepted.  Formulation  of  t(N),  for  the  case  where  f ^  and 

f  are  Gaussian,  leads  to  t(N)  =  T.  [^(^)  ~  / ,(^^3 — 

°  k=l  l}  + 

where  y(k)  is  the  Bayes  estimate  of  y(k)  obtained  recursively  [l,2], 

2  2 

a  is  the  noise  variance,  and  L  is  the  estimated  signal  variance.  The 
required  threshold  can  be  obtained  from  specification  of  error  proba¬ 
bility  (false  alarm  and  miss)  and  p(t  |h^)  p(t  |h^),  where  p  is  the 
probability  distribution  of  t  given  H.,  j=0,  1.  The  errors  corresponding 
to  wrong  decisions  will  then  be  found  from: 


false  alarm  =  J  p(t  jH  )  dt  and 
t  ° 

error  of  second  kind  =  P(t  |Hj)  dt. 

The  above  procedure  has  been  applied  to  a  number  of  pictures 
under  very  low  signal  to  noise  conditions  (peak  to  peak  signal  level/vari¬ 
ance  of  noise  of  about  ^  ).  This  is  far  below  the  level  where  a  human 
observer  can  detect  the  presence  of  signals.  The  preliminary  results  have 
been  very  encouraging.  Attempts  are  being  made  to  use  the  detection 
process  in  conjunction  with  estimation  to  improve  the  image  enhancement 
process  beyond  what  is  achievable  by  linear  estimation. 
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5.2  Edge  Detection  Techniques 

Charles  B.  McGregor  and  Richard  P.  Kruger 

A  fundamental  process  in  pattern  recognition  and  image  process¬ 
ing  is  the  detection  of  edges  in  a  scene.  An  investigation  of  edge  detec¬ 
tion  techniques  appli'id  to  continuous  motion  images  has  been  undertaken. 
The  problem  undei  study  consist.^  of  identifying  and  measuring  the  area 
of  a  brightened  region  on  an  image  frame  and  tracking  this  area  as  it 
changes  from  frame  to  frame.  The  work  performed  to  date  has  con¬ 
sisted  of  an  examination  of  differentiator  and  tracing  types  of  edge  de¬ 
tectors. 

Differentiating  Edge  Detectors.  The  computation  of  edge  elements  by 
differentiation  on  a  discrete  image  takes  the  form  of  simple  additions  and 
substractions  of  pixel  intensities.  The  simplest  such  differencing  oper¬ 
ation  takes  the  form: 


d.  .  =  a.  .  -  a.  ,  . 
i,J  i,J  1-1, J 

Other  forms  which  have  been  used  include  the  gradient  [  1, 2] 

2  2 

d.  .  =  (a.  -a.  .)  +  (a.  .-a.  .) 

i.J  i,J  +  l  i,J  1-1, J  i,J 

and  the  Roberts  gradient 

d.  .  -  |a.  .-a  |  +  |a.  ...-a  .| 

1,  J  1,  J  1+1,  J  +  1  1,  J  +  1  1+1,  J 

Variants  on  these  basic  approaches,  including  the  Laplacian,  have  been 
reported  in  the  literature. 

A  fundamental  problem  occurring  with  simple  differencing  techniques 
is  that  in  real  images  large  differences  may  occur  on  a  pixel  to  pixel  basis 
at  points  both  interior  and  exterior  to  a  desired  edge.  This  causes  dif¬ 
ficulty  if  a  threshold  criterion  is  to  be  used  to  collect  the  set  of  edge  points 
from  the  differentiated  scene.  One  approach  to  this  problem  is  the  use  of 
heavy  smoothing  in  the  computation  of  the  differences,  eg: 
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i.J  1-4, j  1-1. J  1.3  1+3,3 

where  the  superscript  indicates  the  number  of  elements  included  in  each 
of  the  averaging  sums.  A  further  refinement  involves  the  use  of  a  pro¬ 
duct  of  terms: 

,  _  ,1  .  ^2  ,^4  . 

d.  .  -  d  d  d  ... 

1.3 

While  this  product  edge  detector  does  succeed  in  removing  many  unde¬ 
sired  detections  from  the  results,  it  does  so  at  the  price  of  lost  resolu¬ 
tion  and  the  generation  of  a  set  of  sequencies  at  which  the  detector  is 

totally  blind  to  periodic  edges.  The  blind  sequency  effect  can  be  seen 

N 

from  Figure  1.  Consider  the  edge  detector  d.  .  .  As  this  detector  is 

1.  3 

translated  along  the  x-axis,  it  can  be  observed  that  any  waveform  which 
causes  the  integral  under  the  positive  and  negative  cusps  of  the  detector 
to  vary  identically  will  not  contribute  to  the  detector  output.  There  exists 
a  broad  class  of  periodic  waveforms  having  period  equal  to  N/k  (k  an 
integer)  which  possess  this  property.  Thus,  for  example,  the  detector 

g 

d  is  blind  to  waveforms  with  sequencies  of  8,  16,  32,  .  .  .  While  such  an 
effect  is  acceptable  in  situations  in  which  the  data  is  known  to  have  no 
spectral  components  at  the  blind  sequencies,  it  is  totally  unacceptable  for 
cases  in  which  the  spectrum  of  the  data  is  unknown.  While  further  pro- 
ceedures  might  be  developed,  possibly  using  different  values  of  N  with 
suitable  combining  logic,  the  growing  complexity  of  the  detector  becomes 
unattractive  for  an  operation  which  must  be  repeated  on  each  pixel  of  a 


In  the  presence  of  noise,  the  performance  of  the  differentiating 
detector  is  severely  degraded.  Noise  will  increase  the  number  of 
randomly  distributed  points  which  cross  the  derivative  threshold,  causing  a 
rough,  thickened  edge  region  rather  than  the  desired  thin  line.  As  the 
differentiating  detector  does  not  label  interior  points  in  the  desired  set,  edge 
tracer  techniques  as  described  in  the  next  section  are  not  implemented  as 
effectively  as  with  threshold  detection  processing. 

Threshold- Tracing  Edge  Detectors.  The  determination  of  the  edge 
of  a  set  defined  by  threshold  criteria  begins  with  the  comparison  of  each 


pixel  with  a  threshold.  Pixels  below  the  threshold  values  are  set  to  zero 
and  those  exceeding  the  threshold  are  set  to  a  uniform  value.  Thus  the 
data  is  initially  compressed  to  a  one-bit  representation  on  the  intensity 
axis  before  further  processing.  The  threshold  for  this  comparison  may 
be  either  externally  specified  or  computed  by  seeking  an  extremum  of  the 
histogram  of  the  scene.  The  output  from  this  process  will  contain  iso¬ 
lated  and  discontinuous  sets  of  points,  similar  to  the  output  of  a  differ¬ 
entiating  edge  detector.  At  this  point,  continuity  of  the  edge  is  invoked 
and  an  algorithm  which  seeks  connected  points  on  the  frontier  of  the  set 
is  used.  In  topological  terms,  this  algorithm  seeks  points  which  (1)  are 
adjacent  to  previously  determined  perimeter  points  and  (2)  are  the  centers 
of  unit  rr.dius  neighborhoods  containing  points  of  both  the  circumscribed 
set  and  its  complement. 

The  operation  of  this  algorithm  is  best  described  by  ascribing 
its  attributes  to  an  intelligent  bug  which  walks  about  the  image  plane  per¬ 
forming  the  tracing  [3].  Initially,  the  bug  is  either  placed  on  a  known 
boundary  point  or  it  is  caused  to  enter  the  scene  along  a  suitable  horizontal 
or  vertical  line  and  proceeds  until  it  identifies  an  edge  point.  One  typical 
procedure  defines  the  initial  edge  point  as  the  first  interior  point  of  the 
set  which  is  followed  by  a  continuous  line  of  N  points  also  belonging  to 
the  interior.  When  thus  initialized,  the  bug  looks  toward  the  interior  of 
the  set  and  begins  searching  a  unit  radius  path  in  a  clockwise  direction 
until  it  locates  a  point  belonging  to  the  complement  of  the  set.  Having 
found  such  a  point,  the  bug  continues  its  scan  from  that  point  until  it 
encounters  the  first  member  of  the  desired  set.  The  bug  then  moves  to 
the  new  point  and  enters  it  into  a  list  of  perimeter  points.  Next,  it  looks 
back  to  the  previous  point  and  begins  the  search  for  an  exterior  point,  re¬ 
peating  the  above  procedure.  The  bug  repeats  this  procedure,  moving  in 
a  clockwise  direction  around  the  desired  set  until  it  returns  to  the  initial 
point.  Having  traced  a  closed  path,  the  bug  retires.  It  is  noted  that  the 
bug  lives  in  a  square  world  and  considers  the  eight  elements  shown  in 
Fig.  2  as  belonging  to  a  neighborhood  of  unit  "radius".  Since  the  bug 
moves  in  fixed  length  steps,  a  complete  record  of  its  path  consists  of 
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the  coordinates  of  the  initial  point  and  the  series  of  directions  it  moves 
on  each  succeeding  step.  This  highly  compact  representation  of  the 
perimeter  is  known  as  a  chain  code  [4].  This  encoded  data  contains  all 
information  necessary  to  compute  the  enclosed  area  and  the  moments  of 
the  traced  centour  [sj. 

In  the  presence  of  noise,  the  hug  algorithm  performs  admirably 
since  it  only  senses  and  reacts  to  noise  elements  which  lie  within  a 


unit  radius  of  a  perimeter  point.  Further  the  probability  that  the  bug 
will  deviate  from  the  desired  perimeter  by  N  units  equals  the  pro¬ 
bability  that  a  chain  of  N  adjacent  noise  elements  will  occur  with  one  end 
adjacent  to  the  perimeter.  This  probability  varies  with  the  Nth  power 
of  the  probability  of  a  single  noise  element.  The  only  restriction  on  the 
tracing  of  closely  adjacent  set  perimenters  is  the  probability  the  noise 
will  bridge  between  the  perimenters,  causing  two  sets  to  appear  as  one. 

An  Experimental  Edge  Detection  Study.  The  performance  of  the 
threshold-tracing  edge  detector  has  been  evaluated  using  radiographic 
films  of  the  human  heart.  The  cardiac  catheterization  procedure  used  in 
the  evaluation  of  severe  cases  of  heart  disease  involves  injection  of  radio¬ 
opaque  dye  into  the  left  ventricle  (the  output  pumping  chamber)  of  the  heart 
while  motion  pictures  are  made  of  the  flouroscopic  image  of  the  heart. 
From  these  films,  the  opaque  image  of  the  ventricle  can  be  traced  to  show 
its  contraction  as  a  function  of  time  while  the  heart  beats.  From  the  cross 


sectional  area  shown,  the  volume  may  be  estimated  as  a  function  of  time 


and  the  volume  of  blood  expelled  by  the  heart  per  beat  can  be  computed  for 


diagnostic  use.  At  a  motion  picture  rate  of  thirty  frames  per  second,  a 
manual  processing  of  these  fifty  images  requires  ten  to  twelve  hours  of  a 
skilled  technician's  time.  Thus  a  need  for  computer  data  processing  arises. 

The  experimental  data  consists  of  strips  of  35  millimeter  film  which 
have  been  scanned  in  a  flying  spot  scanner.  A  view  of  the  heart  with  no  dye 
injected  is  available  to  allow  subtraction  of  background  opacity  from  ihe 
frame  to  be  processed.  A  frame  of  film  showing  a  ball  of  known  diame'er 
in  the  position  later  occupied  by  the  heart  is  available  for  calibration  pur¬ 
poses. 
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Initial  processing  consists  of  subtraction  of  the  background 
picture  from  the  first  image  to  be  processed  and  computation  of  a 
gray  level  histogram  for  the  resulting  p'^-iure.  A  bimodal  histogram 
is  assumed  and  a  curve  is  fitt.d  to  the  valley  'if  the  histogram.  The 
central  minimum  of  this  curve  is  used  as  the  threshold  for  edge  de¬ 
tection.  After  extensive  investigation,  a  version  of  the  bug  which  does 
not  require  prior  thresholding  of  the  data  has  evolved.  This  bug  ac¬ 
cepts  threshold  levels  from  the  calling  program  and  makes  the  nec¬ 
essary  threshold  comparisons  as  it  traces  its  path  through  a  field 

of  unthresholded  data.  This  procedure  reduces  the  threshold  com- 

2 

parisons  from  N  to  approximately  16N  operations  and  allows  the  data 
to  be  viewed  with  the  threshold  tracing  superimpo.sed.  A  typical  trac¬ 
ing  is  accomplished  in  0.  8  seconds.  Examples  of  raw  data  and  the  data 
with  an  edge  tracing  superimposed  are  shown  in  figures  3  and  4. 

Points  of  Departure.  The  next  task  in  applying  the  contour  tracing  to 
frame  to  frame  tracking  involves  reducing  the  large  amount  of  data 
present  in  a  chain-encoded  contour  to  a  few  descriptive  paramenters. 

An  initial  approach  will  use  a  Fourier  expansion  on  the  closed  boundary. 

The  minimum  number  of  Fourier  coefficients  needed  to  approximate  the 
boundary  within  an  area  erro'  limit  will  be  established.  A  parallel 
investigation  is  underway  which  may  result  in  the  integration  and  chain 
code  processing  portions  of  this  project  being  used  with  manual  Graf- 
pen  input  of  data  to  a  time-sharing  computer  system  for  interim  man- 
machine  data  processing. 
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6.  Image  Processing  Support  Projects 

The  following  sections  describe  the  status  of  the  U.  S.  C. 
image  processing  research  facilities  with  particular  emphasis  on 
thp  pr'^gress  iov/ard  rnakiiig  the  image  processing  hardware  and 
software  facilities  available  over  the  ARPANET.  A  report  is  also 
given  on  the  initial  steps  towards  the  establishment  of  a  standard 
set  of  original  images  for  the  image  processing  research  commu¬ 
nity. 

6.1  U.  S.  C.  Image  Processing  Hardware  and  Software  Facilities 
Richard  P.  Kruger 

Image  Processing  Hardware  Facilities.  Figure  1  is  a  gen¬ 
eral  block  diagram  depicting'  the  hardware  facilities  utilized  by  the 
U.  S.  C.  Image  Processing  Laboratory  (IPL)  for  various  image  acqui¬ 
sition,  computation,  and  display  operations.  The  IPL  facility  con¬ 
sists  of  a  Hewlett-Packard  2100  which  is  used  to  command  and  con¬ 
trol  the  image  acquisition  and  display  devices,  and  an  Interdata  IV 
which  is  used  for  simulation  of  image  processing  systems. 

At  present  a  conventional  flying  spot  scanner  is  used  for  image 
digitization  and  hard  copy  display  of  color  and  monochrome  images  on 
photographic  prints.  The  device  can  raster  scan  transparencies  up  to 
70  millimeters  square  with  resolutions  up  to  1024  by  1024  pixels.  The 
scan  time  varies,  bur.  is  usually  less  than  1  minute  per  frame.  Normal 
analog  to  digital  conversion  is  eight  bits.  Display  is  usually  on  Polaroid 
or  35  mm  film.  Color  digitization  and  display  is  accomplished  serially 
under  computer  control  using  a  color  filter  wheel. 

The  Muirhead  color  facsimile  device  is  a  two  unit  drum  trans¬ 
mitter  and  receiver.  Its  capabilities  include  100  line  per  inch  resolu¬ 
tion  with  simultaneous  ten  bit  per  color  quantization  as  well  as  display 
of  the  three  primary  colors  on  an  8  by  10  inch  print.  The  scan  time  for 
this  device  is  12  minutes.  While  the  maximum  scan  or  display  size  is 
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Figure6.  1-1^  USC  Image  Processing  Laboratory  Facilities 


fixed  on  an  8  by  10  inch  field,  options  exist  for  scan  or  display  of  any 
2,  5  or  5.  0  inch  square  portion  of  a  larger  field  at  resolutions  of  256  by 
256  or  512  by  512  pixels. 

The  real  time  display  of  monochrome  and  color  digital  images 
is  accomplished  using  an  Aerojet  General  display  device.  This  device 
utilizes  a  standard  shadow  mask  CRT  with  576  horizontal  and  525  ver¬ 
tical  lines  of  resolution.  A  digital  disk  is  included  with  the  device  which 
makes  possible  a  refresh  rate  of  60  fields  per  second  at  64  quantization 
levels  for  each  of  the  red,  green,  and  blue  primaries. 

An  Optronics  flat  bed  scanning  microdensitometer  is  scheduled 
for  delivery  by  1  January  1973.  A  specification  of  the  unit  is  enclosed 
as  an  appendix  to  this  report.  The  Optronics  unit  will  scan  and  display 
16,  35,  and  70  millimeter  color  or  monochrome  roll  film  on  a  registered 
transport  under  computer  control  as  well  as  provide  prints  or  transparen¬ 
cies  up  to  10  inches  square.  Specifications  call  for  a  minimum  aperture 
size  of  2  microns  square,  12  bits  or  precision,  a  scanning  velocity  of  8 
inches  per  second  in  the  direction  of  travel  and  a  0  to  4  specular  density 
range.  Color  digitization  and  display  is  sequential  using  a  color  wheel  un¬ 
der  program  control. 

The  IPL  computers  and  peripherals  are  hardware  interfaced  to 
the  Engineering  Computer  Laboratory  (ECL)  IBM  360/44.  Images  may  be 
transferred  from  the  disk  on  the  HP  2100  to  a  disk  on  this  machine,  pro¬ 
cessed  and  then  transferred  back  to  the  HP  2100.  Work  is  continuing  on 
an  operating  system  that  will  provide  remote  job  entry  between  the  HP  2100 
and  the  IBM  360/44  for  image  processing  tasks.  The  IBM  360/44  is  opti¬ 
cally  linked  to  the  University  Computer  Center  (UCC)  IBM  370/155  for  re¬ 
mote  job  entry  between  the  machines.  The  IBM  360/44  also  acts  as  the 
U.  S.  C.  host  computer  for  the  ARPANET  using  a  TIP  interface. 

Image  Processing  Software  Capability.  Image  processing  tasks 
on  the  IBM  370/155  are  performed  with  conventional  FORTRAN  based  soft¬ 
ware  and  a  processing  software  package  called  VICAR.  The  VICAR  sys- 
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tem  is  also  being  implemented  on  the  IBM  360/44.  The  VICAR  (Video 
Image  Commanication  eind  Retrievel)  system  was  acquired  from  the  Jet 
Propulsion  Laboratory,  and  a  inodifiad  form  has  been  operational  at  Q.  S.  C. 
for  the  past  year. 

The  purpose  of  VICAR  is  to  facilitate  the  acquisition,  digital  pro¬ 
cessing,  eind  recording  of  image  data  on  a  production  basis  by  scientific 
personnel  not  familiar  with  systems  programming,  applications  program¬ 
ming,  and  problem  programming.  The  system  programmer  requires  know¬ 
ledge  of  applications  and  problem  progranriming  only  in  a  global  sense.  The 
applications  programmer  is  capable  only  of  using  existing  problem  programs 
and  requires  no  knowledge  of  either  system  or  problem  programming.  The 
results  of  processing  an  image  with  problem  programs,  GRID,  MAPGRID 
and  CONCAT  are  shown  in  Figure  2a,  b  and  c.  The  problem  programmer 
must  have  knowledge  of  application  program  protocol  but  need  not  have 
system  programming  knowledge.  From  a  system  viewpoint,  the  opera¬ 
tion  of  VICAR  is  very  similar  to  that  of  a  standard  utility  program.  The 
application  programmer  submits  a  card  deck  containing  a  limited  number 
of  VICAR  control  statements  to  define  the  processing  tasks.  These  state¬ 
ments  include  READ  and  WRITE  tape  instructions  in  7  or  9  track,  6  or  8 
bit  formats,  disk  data  set  RESERVE  cards  and  program  EXECUTE  cards. 

All  tape  and  disk  labels  and  job  control  statements  are  generated  so  as  to 
be  transparent  to  the  user.  Up  to  14  simultaneous  I/O  requests  can  be  made. 
User  labels  as  well  as  history  and  system  labels,  are  routinely  placed  at 
prescribed  locations  preceding  the  image  data  on  a  tape  or  disk  file.  This 
information  can  be  retrieved  and  used  for  picture  annotation  purposes  such 
as  shown  in  Figure  2d.  This  is  one  display  option  of  the  problem  program 
MASK.  Labels  can  be  deleted  in  total  or  in  part  at  the  user's  option.  The 
problem  programmer  unlike  the  applications  programmer  can  add  programs 
to  the  program  library  which  can  be  executed  by  himself  or  by  an  appli¬ 
cations  programmer.  A  problem  program  can  be  written  either  in  Fortran 
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(a)  Output  of  GRID 


(b)  Output  of  MAPGRID 


(c)  Output  of  CONG  AT 


(d)  Output  of  MASK 


(f)  Imapc  processed  by  Gaus 
sian  isotropic  high  frequency 
phasis  filter  followed  by  histo 
gram  equalization. 


Examples  of  Processed  Images 


or  Assembly  Language  and  differs  from  a  normal  subroutine  only  in 
its  I/O  structure.  In  addition  to  i.ornrial  Fortran  I/O  which  is  option¬ 
al,  several  VICAK.  I/O  staienients  must  be  included.  Most  existing 
conventional  subroutines  can  be  easily  converted  to  VICAR  problem  pro¬ 
grams  through  a  modification  of  its  I/O  structure.  At  present  VICAR 
contains  over  30  application  programs  which  perform  such  varied  tasks 
as  Fourier  transformations,  picture  annotation,  line  plotting  and  filter 
design.  An  updated  version  of  VICAR  with  block  I/O  transfer  has  been 
installed  at  UCC.  This  modification  will  decrease  I/O  time  and  signi¬ 
ficantly  reduce  I/O  costs  involved  in  image  processing  on  that  machine. 

APPENDIX 

Specifications  for  Proposed  Microdensitometer 
Film  Scan  and  Playback  System 
for 

University  of  Southern  California  Image  Processing  Laboratory 

A.  Digitizer/ Display  Specifications 

1.  The  ability  to  scan  and  playback  sprocketed  roll  film  transparen¬ 
cies  of  16,  35,  and  70  millimeter  sizes  using  an  index  registered 
automatic  film  transport  which  may  be  advanced  or  reversed  manu¬ 
ally  and  under  program  control  an  arbitrary  number  of  frames. 

The  transport  device  must  provide  frame  to  frame  registration 
error  of  not  greater  than  +  12  microns  perpendicular  and  parallel 
to  direction  of  film  travel. 

2.  The  ability  to  directly  digitize  and  playback  sheet  film  photographic 
transparencies  of  opacities  (prints)  up  to  a  maximum  size  of  25.6 
cm  square. 

3.  Color  separation  by  means  of  a  3  primary  color  wheel  under  manual 
and  program  control  to  digitize  or  playback  color  images  in  a  se¬ 
quential  mode. 

4.  The  capacity  to  digitize  and  playback  images  in  a  range  from  0  to 
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4  specular  sensity  units  (0  -  100%  transmittance)  with  at 
least  12  bit  A/D  and  10  bit  D/A  conversion  with  at  least  8 
bits  of  significance  minimum  in  any  application.  |Of  course 
film  has®  fog  and  maximum  density  sensiti\'ity  wi'\  affect 
particular  applications.) 

5.  Linearity  of  +  .02  specular  density  units  (+  .5  percent  trans¬ 
mittance)  over  the  above  stated  range  and  half  the  above  varia¬ 
tion  over  any  2  density  unit  window  within  the  0-3  density  unit 
range. 

6.  Stability  of  less  than  ^  .02  density  unit  drift  .  5  percent  trans¬ 
mittance)  after  one  hour  of  warm  up  to  be  sustained  for  at  least 
10  hours. 

7.  Repeatability  of  ^  .01  density  units  or  less  after  one  hour  of 
warm  up. 

8.  Optical  system  depth  of  focus  to  exceed  250  microns  at  all  aper¬ 
ture  sizes  or  a  means  whereby  the  transported  roll  film  or  mounted 
sheet  film  curvature  will  not  exceed  the  depth  of  focus  at  the  film 
plane  for  any  selected  aperture. 

9.  Square  apertures  capable  of  256,  512,  and  1024  line  resolutions  for 
a  square  frame  scan  or  replay  on  all  types  roll  film  as  well  as  for 
sheet  film  sized  of  4,  5,  8,  and  10  inches  square  with  no  aperture 
overlap.  (It  is  expected  that  standard  apertures  available  with  the 
device  will  also  be  included  and  that  other  apertures  may  be  easily 
added  in  the  future, ) 

10.  X  and  Y  axis  spatial  repeatability  of  less  than  or  equal  to  1 
micron. 

11.  X  and  Y  axis  straightness  of  travel  variation  less  than  or  equal  to 
^  5  microns. 

12.  X  and  Y  axis  orthogonality  of  ^  .  5  seconds  of  arc. 

13.  Position  accuracy  of  1  micron  RMS  in  X  and  Y  direction  in  25.6  cm 
of  travel.  (6  micron  peak  to  peak) 
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14.  Digitization  or  replay  speed  variable  from  approximately  5  to 
20  cm/ second  minimal. 

15.  Peak  signal -to-noise  ratio  exceeding  50  db  in  large  aperture 
clear  field  operation  and  to  exceed  40  db  under  all  opt-xating 
conditions. 

16.  Photomultiplier  surface  material  must  be  capable  of  extended 
red  wavelength  response  to  provide  minimum  signal-to-noise 
ratio  requirement  for  color  image  scanning.  It  is  suggested 
that  S-20,  S-25,  or  Gallium  Arsenide  photosurfaces  with  a  thermal 
electric  PMT  cooler  be  investigated  to  meet  this  specification. 

B.  Other  Specifications 

1.  A  viewing  screen  for  accurate  manual  film  positioning  as  well  as 
a  means  for  viewing  small  portions  of  the  film  plane. 

2.  Switchable  logarithmic  converter. 

3.  Direct  digital  display  of  film  density  information. 

4.  Console  display  with  convenient  switches  and  indicators  such  as 
end  of  scan  light,  etc. 

5.  Means  for  accurate  initial  and  interval  calibration  of  instrument  to 
ensure  continuing  densitometer  accuracy. 

6.  The  proposed  system  will  be  interfaced  to  a  Hewlett-Packard  2100 
with  16  bits  in  and  out,  flag  in  and  out  configuration,  A  typical 
I/O  card  for  this  purpose  is  the  microcircuit  interface  card  No. 
12566.  Cooperation  with  our  computer  staff  and  Hewlett-Packard 
will  be  available  to  aid  in  software  development  and  consultation. 

7.  Switching  circuitry  to  allow  bit  packing  into  a  density  or  transmit¬ 
tance  range  existing  within  a  particular  transparency  or  print  so 
that  the  potential  exists  for  full  register  length  usage. 

8.  A  complete  maintenance  manual  so  that  our  laboratory  staff  can 
perform  all  repairs  and  periodic  adjustments. 

9.  Scan  and  replay  capability  in  a  normal  room  light  environment. 
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6.2  use/ ARPANET  Image  Processing  System 
Richard  P.  Kruger  and  William  K.  Pratt 

The  University  of  Southern  California  image  processing  hard¬ 
ware  and  •jort-wii.rc  facilities  previo isiy  discussed  will  soon  be  available 
for  use  over  the  ARPA  computer  network.  A  remote  site  user  may: 

(a)  read  an  image  from  an  image  file  and  transfer  the 
image  to  the  remote  site; 

(b)  store  an  image,  transferred  from  the  remote  site, 
on  an  image  file; 

(c)  display  and  make  a  hard  copy  of  an  image  stored  on 
an  image  file; 

(d)  digitize  a  hard  copy  image  and  store  on  an  image  file; 

(e)  execute  VICAR  language  image  processing  program. 

The  use  of  the  USC/ARPANET  image  processing  system  by  a  remote  site 
user  requires,  as  a  minimal  configuration,  a  computer  terminal  connected 
to  a  TIP.  With  this  configuration  image  processing  jobs  can  be  initiated 
from  the  remote  site  with  operations  performed  on  previously  stored  images 
and  hard  copy  image  results  returned  to  the  user  by  mail.  The  next  level 
of  capability  entails  the  addition  of  an  image  recording  device  at  the  re¬ 
mote  site  for  display  of  processed  images.  A  complete  remote  site  sys¬ 
tem  would  also  include  an  image  scanner  and  digitizer  for  supplying  input 
images. 

Figure  1  contains  a  block  diagram  describing  the  information  flow 
in  the  USC/ARPANET  image  processing  system.  In  operation  a  user  at  a 
remote  site  establishes  a  connection  with  the  USC  network  control  program 
and  notifies  the  NCP  that  an  image  processing  job  is  to  be  executed.  A 
data  transfer  program  receives  a  request  for  an  image  processing  job  from 
the  NCP  and  then  acts  as  a  monitor  for  control  of  the  job  instructions  and 
data.  The  remote  user's  job  parameters,  consisting  of  instructions  defin¬ 
ing  the  image  processing  operations  to  be  performed  and  image  data  loca- 
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Figure  6.2-1  USC/ARPANET  IMAGE  PROCESSING  SYSTEM 


tions,  are  stored  in  a  job  file  preparatory  to  execution.  If  the  user  sub¬ 
mits  a  VICAR  program,  the  source  program  is  stored  on  a  mass  iniage 
file.  When  the  remote  user’s  job  is  ready  for  execution  the  data  transfer 
program  transfers  the  appropriate  image  data  and  instruction  to  the  pro¬ 
cessors.  For  example,  if  a  ^^JCAR  job  if  to  be  iie.rformed  the  data  tsaus- 
fer  program  transfers  the  remote  user's  VICAR  language  source  program 
and  the  input  images  to  be  processed  from  the  image  files  to  the  VICAR 
language  processor  on  the  IBM  370/155.  After  the  VICAR  processing  is 
completed,  the  output  image  is  returned  to  the  image  file,  and  subse¬ 
quently  returned  to  the  remote  site  through  the  TIP. 

The  majority  of  messages  transmitted  over  the  network  are  of 
short  length  and  can  be  placed  in  a  single  packet.  However,  bulkfiles 
such  as  digital  images  usually  exceed  the  packet  length,  and  therefore, 
must  be  partitioned  for  transmission.  The  partitioning  is  automatically 
handled  by  the  host  and  user  network  control  programs,  and  is  not  of  con¬ 
cern  to  the  user.  For  digital  image  processing  operations,  images  are 
usually  coded  with  eight  bits  (one  byte)  ner  picture  element.  Further¬ 
more,  images  are  often  restricted  in  size  to  a  binary  integer,  e.  g.  (256 
by  256,  512  by  512,  etc.)  in  order  to  take  advantage  of  fast  computational 
algorithms  which  impose  this  restriction.  In  the  interest  of  standardiza¬ 
tion,  the  USC/ARPANET  image  processing  system  has  been  designed  to 
handle  image  files  with  a  record  length  of  256N  bytes  per  record  where 
N  is  an  integer.  The  number  of  records  in  a  file  can  be  specified  by  the 
user,  and  is  not  restricted  to  a  binary  integer.  At  present  the  fastest 
transmission  speed  over  the  ARPANET  is  501  bs.  Considering  the  trans¬ 
mission  control  data  that  automatically  is  added  to  each  message  packet 
and  transmission  delays  encountered  at  each  node,  the  average  user  data 
rate  is  likely  to  be  about  25kbs.  At  this  rate  the  time  required  to  trans¬ 
mit  a  256  by  256  byte  image  is  about  21  sec.  A  1024  by  1024  image  re¬ 
quires  about  5.6  min,  for  transmission.  Plans  are  underway  to  provide 
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a  higher  speed  transmission  service  to  the  network.  Also,  image  band¬ 
width  reduction  techniques  are  being  considered  to  decrease  the  transmis¬ 
sion  time. 


6.  3  Development  of  Original  Images  for  Image  Processing 
William  K.  Pratt 

A  master  set  of  three  monochrome  and  two  natural  color  images 
has  been  assembled  as  a  provisional  standard  for  image  coding  studies. 

The  following  sections  describe  the  original  master  set.  It  is  antici¬ 
pated  that  future  master  sets  will  contain  a  wider  variety  of  imagery, 
useful  for  image  coding,  restoration,  enhancement,  etc.  studies,  and 
that  the  subject  pictures  will  be  better  calibrated,  analyzed,  and  documented. 

General  Picture  Description.  All  test  pictures  contain  256  by  256 
pixel  points.  The  luminance  (intensity)  at  each  pixel  point  has  been  li¬ 
nearly  quantized  into  256  levels  and  binary  coded  into  the  range  (0,  255). 

The  test  pictures  are  listed  below: 


File  # 


Name 


1  girl  -  red  tristimulus  value 

2  girl  -  green  tristimulus  value 

3  girl  -  blue  tristimulus  value 

4  couple  -  red  tri stimulus  value 

5  couple  -  green  tristimulus  value 

6  couple  -  blue  tristimulus  value 

7  girl  -  monochrome  luminance 

8  couple  -  monochrome  luminance 

9  moonscape  -  monochrome  luminance 
The  monochrome  luminance  picture  of  the  girl  and  couple  have  been  ob¬ 
tained  from  the  natural  color  representation  through  the  coordinate  con¬ 
version 

=  .299R^  +  .  587G^  +  .  114B^ 
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where  =  image  luminance  in  integer  form  (0  <  <  255) 

R^,  G^,  =  red,  green,  blue  image  tristimulus 

value  in  integer  form  (0  <  <  255), 

255).  (0  <  B^<  255). 

The  moonscape  image  is  a  picture  of  the  moon  taken  from  a 
spacecraft  and  subsequently  processed  by  the  Jet  Propulsion  Laboratory 
to  enhance  its  subjective  qulaity.  No  information  is  readily  available 
as  to  characteristics  of  the  image  sensor  or  the  post  processing  algo¬ 
rithms. 

The  color  images  have  been  scanned  by  Kodak  from  color  slides 
distributed  by  the  Society  of  Motion  Picture  and  Television  Engineers 
(SMPTE).  The  digital  tapes  and  supporting  documentation  were  supplied 
to  use  by  Dr.  Oleh  Tretiak  of  M. I.  T.  In  the  scanning  process,  the 
scanner  output  (nominally  photographic  density)  was  linearly  quantized 
to  256  levels.  System  calibration  curves  were  then  generated  to  relate 
a  scanner  digital  number  to  true  density  and  true  transmittance.  At 
use  the  calibration  curves  were  used  to  convert  the  scanner  digital  num¬ 
bers  to  a  set  of  digital  numbers  linearly  proportional  to  the  transmittance 
of  the  test  slides  for  each  tristimulus  value.  The  resulting  transmission 
scale  between  0.  and  0.  6  has  been  linearly  quantized  to  256  levels.  The 
SPMTE  slides  were  scanned  over  an  equi- space  raster  of  520  lines  and 
780  pixels  per  line.  The  couple  picture  of  the  master  set  was  obtained 
from  the  scan  of  SMPTE  set  #2  by  extracting  the  center  512  by  512  pixels 
and  spatially  averaging  the  corrected  transmittance  in  groups  of  four 
pixels  to  obtain  a  256  by  256  pixel  image.  The  girl  picture  of  the  master 
set  was  obtained  by  extracting  a  256  by  256  pixel  area  of  the  same  digi¬ 
tized  SMPTE  slide  #3. 

Statistical  measurements  have  been  performed  on  the  master 
set  ima'fes.  These  measurements  include: 
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mean,  variance,  standard  deviation 
row  and  column  covariance 
grey  scale  histogram  (64  levels) 
first  order  entropy 

second  order  row  and  column  entropy 

For  storage  of  images  on  magnetic  tape,  the  tape  format  is 

2400  ft.  reel 
9  track 
800  bpi 

one  pixel  per  byte 
256  bytes  per  record 
256  records  per  file 
no  tape  label 

record  gap  after  each  record  (IBM  compatable) 
file  mark  after  each  file  (IBM  compatable). 
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7.  New  Resejarch  Projects 


The  following  are  summary  descriptions  of  new  research  studies 
that  will  be  initiated  during  the  next  six  months  in  addition  to  the  continu¬ 
ation  of  topics  presented  earlier. 

7.  1  Interframcj  Coding  of  Television  Signals 
Ali  Habibi 

The  contour  tracing  algorithm  has  been  employed  as  a  technique 
of  encoding  two  dimensional  digital  data  by  tracing  the  contours  of  equal 
gray  levels  and  transmitting  the  addressing  information  that  enables  the 
receiver  to  reconstruct  these  contours.  This  technique  relies  upon 
the  spatial  correlation  to  reduce  the  data  redundancy.  It  is  proposed  to  use 
the  contour  tracing  algorithm  to  encode  interframe  differential  pictures. 
These  pictures  are  created  by  taking  the  difference  of  two  successive 
frames  in  the  television  signal  and  indicating  the  change  along  the  tem¬ 
poral  axis.  The  proposed  technique  reduces  the  redundancy  in  both 
sptaial  and  temporal  domains,  therefore  achieving  a  further  reduction 
in  the  redundancy  than  is  possible  with  either  one  of  the  techniques.  Its 
potential  use  includes  coding  both  monochrome  and  color  video  signals. 

Preliminary  results  indicate  that  good  coding  capabilities  may  be 
achieved  by  combining  unitary  transformations  with  a  DPCM  encoder. 

Using  this  technique  for  intraframe  coding,  first  a  one  dimensional  uni¬ 
tary  transformation  is  used  to  obtain  the  transform  of  each  line  of  a 
still  picture.  This  gives  a  signal  in  the  transform  domain  which  has 
a  strong  correlation  in  the  vertical  direction,  then  this  signal  is  coded  by 
a  DPCM  coder  taking  advantage  of  this  correlation. 

It  is  planned  to  extend  this  technique  to  interframe  coding  by  taking 
the  two-dimensional  unitary  transform  of  each  frame  and  then  applying 
the  DPCM  encoder  to  code  the  transformed  data  in  the  temporal  direction. 
Such  a  system  will  explore  the  correlation  of  moving  pictures  in  spatial 
as  well  as  temporal  directions  and  allow  a  parallel  operation  that  can  be 
used  to  implement  real-time  interframe  coders. 
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7.  2  Space  Variant  Point  Spread  Functions 
Harry  C.  Andrews 

Due  to  the  increase  in  connputing  technology,  it  is  becoming 
feasible  to  consider  the  task  of  computationally  implementing  image 
restoration  systems  degraded  by  space  variant  point  spread  functions 
(SVPSF).  Five  specific  areas  of  investigation  under  this  subject  are: 

1.  Object  dependent  SVPSF' s 

2.  Eigenfunction  Analysis 

3.  Degree  of  Freedom  and  SVPSF's 

4.  Gaussian  Kernel  SVPSF's 

5.  Atmospheric  degradating  SVPSF's 

Each  one  of  these  areas  is  being  investigated  with  respect  to  analysis 
and  computational  implementation.  The  object  dependent  SVPSF  is 
modeled  after  the  equation 


g(x,y) 


f(§,  T])  h(x,  y,  §,  Tl,  f)  d§dTi 


An  approach  to  the  solution  of  this  equation  is  a  Taylor  Series  expan¬ 
sion  of  h(.)  about  f(.).  The  Eigenfunction  analysis  utilizes  the  equation 


=  ff  ♦„J5,n)h(5.r,,x,y)  d5dn 

-  00 

to  develop  a  system  of  Eigenfunctions  (if  they  exist)  for  expansion  of  the 
image  function. 

The  degrees  of  freedom  analysis  results  in  a  four-dimensional 
singular  valued  decomposition  to  desc  '^e  h(x,  y,  Tj)  more  efficiently 
for  computational  inversion.  The  Gaussian  kernel  SVPSF  analysis 
attempts  to  invert  a  point  spread  function  given  by 

h(p.  q)=  A(p)  exp  -  {  (p  -  q)^K  ^(p)(p  -  q)  ] 

where 

£  =  (x,  y) 

a. "  (?.  Tl) 

Various  simplications  are  introduced  for  computational  simplicity. 
The  atmospheric  degradation  SVPSF  is  really  a  phase  sensitive  SVPSF 


-116- 


system  in  which  the  image  and  point  spread  function  phases  are  particularly 
significant.  Phase  unwrapping  and  phase  filtering  will  probably  be  re¬ 
quired  for  such  restoration  processes. 

7.3  Fast  Recursive  Restoration  for  Two  Dimensional  Images 
Anil  K.  Jain 

The  standard  methods  used  for  recursive  image  restoration  have 
the  following  difficulties: 

a)  The  filters  are  not  isotropic,  i.e.  ,  the  improvement  in 
the  S/N  ratio  depends  considerably  on  whether  the  image  is  scanned 
horizontally  or  vertically; 

b)  The  computation  time  for  large  images  is  considerable; 

c)  Vector  scanning  gives  a  better  S/N  ratio  but  increases  dimen¬ 
sionality  which  leads  to  large  storage  requirements. 

A  second  order  vector  Markov  model  for  two-dimensional  images 
has  been  developed  and  it  is  believed  (based  on  some  initial  results)  that 
an  isotropic  filter  with  scalar  filter  equations  (for  vector  scanning)  cen 
be  developed  leading  to  a  fast  recursive  restoration  algorithm  C  1  ]. 
Reference 

1.  E.  Angel  and  A.K.  Jain,  'A  Nearest  Neighbors  Approach  to  Multi¬ 
dimensional  Filtering,”  Conference  on  Decision  and  Control,  New 
Orleans,  December,  197Z;  also  to  ppear  in  IEEE  Transactions  Aut. 
Control,  April,  1973. 

7.4  Pseudocolor  Mapings  Via  Color  Distance 
Anil  K.  Jain 

This  project  under  consideration  will  deal  with  the  determination 
of  hue  and  saturation  distribution  as  functions  of  illuminance  (or  intensity) 
of  a  black  and  white  image  so  that  a  suitably  chosen  criterion  in  terms  of 
color  distance  [  1  ]  is  satisfied. 

Reference 

1.  A.K.  Jain,  "Color  Distance  and  Geodesics  in  Color-3  Space”,  to 
appear  in  the  Journal  of  the  Optical  Society  of  America. 
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7.  5  Image  Registration 

William  K.  Pratt  and  Alexander  A.  Sawchuk 

A  common  image  processing  problem  is  that  of  spatially  reg¬ 
istering  a  pair  of  images  from  the  same  scene.  For  example,  one 
image  may  be  from  a  visible  wavelength  image  sensor  and  the  other 
from  an  infrared  sensor.  The  general  problem  is  to  spatially  reg¬ 
ister  the  images  and  compensate  for  translation,  rotation,  and  mag¬ 
nification  differences  in  the  presence  of  sensor  noise  and  sensor  non- 
linearities. 

A  preliminary  approach  taken  to  the  problem  has  been  the 
development  of  an  N-dimensional  correlation  function  expressing  the 
N  degrees  of  freedom  of  misregistration.  This  technique  has  been 
implemented  by  software  for  translation  errors  with  good  results. 

The  unique  feature  of  the  program  is  its  utilization  of  the  spatial 
correlation  of  each  image  plane. 

Further  studies  will  expand  the  correlation  technique  for  mag¬ 
nification  and  rotation.  Also,  consideration  will  be  given  to  sequential 
search  techniques. 
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